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A study of the dynamics of a tunneling particle in a driven bistable potential which is moderately- 
to-strongly coupled to a bath is presented. Upon restricting the system dynamics to the Hilbert 
space spanned by the M lowest energy eigenstates of the bare static potential, a set of coupled 
non-Markovian master equations for the diagonal elements of the reduced density matrix, within 
the discrete variale representation, is derived. The resulting dynamics is in good agreement with 
predictions of ab-initio real-time path integral simulations. Numerous results, analytical as well as 
numerical, for the quantum relaxation rate and for the asymptotic populations are presented. Our 
method is particularly convenient to investigate the case of shallow, time-dependent potential barri- 
ers and moderate-to-strong damping, where both a semi-classical and a Redfield-type approach are 
inappropriate. 
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I. INTRODUCTION 



The sensitivity of tunneling to the influence of the environment has been in the focus of intense research over the 
last years ■ A popular model for the investigation of tunneling processes is a double- well potential with an energy 
barrier that separates two energetically degenerate minima. In an idealized system, the barrier can be coherently 
traversed by a quantum mechanical particle {coherent tunneling). A real physical system, however, experiences the 
influence of the surrounding "outer world" . This coupling disturbs the coherent tunneling process and it constitutes 
the origin of decoherence and dissipation in the quantum system. To model the dissipative influence, the environment 
is commonly described as an ensemble of harmonic oscillators (heat bath, reservoir) being at thermal equilibrium 
at temperature T. A bilinear coupling between the quantum system and the bath mimics phenomenologically the 
interaction of the system with the "rest of the world" . By this method, the quantum mechanical analogue of the 
generalized Langevin equation can be derived. 

The spectrum of the uncoupled symmetric bistable potential consists of a ladder of doublets being pairs of ener- 
getically nearly degenerate energy eigenstates. The degeneracy is lifted by the tunneling splittings within the single 
doublets. The doublets themselves are separated by large interdoublet energy gaps which are of the order of the 
related characteristic system frequency scale, the latter are generally orders of magnitude larger than the tunneling 
sphttings. 

By now, two different situations have been in the center of detailed investigations on the dissipative tunneling 
dynamics in a bistable potential: (i) On the one hand, one considers the regime of low temperatures, i.e., k^T is 
of the order of the energy splitting of the lowest tunneling doublet. A common approach to simplify the spatially 
continuous dynamics consists then in restricting the problem to the two lowest energy eigenstates, being the solely 
significantly thermally populated in this deep quantum regime. Coupling the two-level system to a bosonic bath of 
harmonic oscillators leads to the prominent spin-boson problem 0]^,^. (ii) On the other hand, the starting point 
is classical rate theory. Semi- classical tunneling corrections to the relaxation rate are calculated by use of various 
instanton techniques Q. This formalism is applicable when the quantized energy levels lie very dense below the 
barrier, i.e., in cases when the energy barrier is large compared to the characteristic level splitting of the quantum 
system. Moreover, a local equilibrium is required, restricting this approach only to time-mdependent systems. By 
complex-time path integral techniques, the free energy is calculated in a semiclassical steepest-descent method. This 
leads to the dissipative bounce solution which in turn determines the semiclassical decay rate. 

Modern experimental developments have paved the way to study the influence of time-dependent external driving 
forces like a laser beam or an rf-field. Such time-dependent driving fields have most interesting implications for 
quantum systems like, for instance, the effect of coherent destruction of tunneling the effect of quantum stochastic 
resonance [|^16|, or the occurrence of quantum steps in hysteresis loops |p^ , ^ , to name but a few (for recent reviews, 
see Jl|,|l|,|]J). Such driving fields may also be used to control and reduce decoherence in open quantum systems 

Mm 

The present work deals with tunneling processes in a time-dependent bistable potential in a temperature regime 
where the two-level approximation {spin-boson regime) is invalid. Likewise, the (possibly) strong time-dependent 
external fields prevents the use of semiclassical methods. Our analysis therefore bridges those two well established 
limiting regimes in quantum rate theory. 

With this objective in mind, we release the restriction of the bistable potential to its two lowest energy eigenstates 
and extend the model to include more energy eigenstates which are populated at higher temperatures. This implies 
an interesting consequence: Since the energy splittings of the higher doublets are larger, tunneling becomes more 
favorable via the higher doublets. However, for the temperature being too large, tunneling is again hampered due to 
the decoherent infiuence of the environment. This interplay among tunneling, vibrational relaxation (i.e., transitions 
between the doublets) and thermal effects leads to a rich and complex dynamics. 

The specific problem we tackle is the following: Let us consider a quantum particle which is initially localized in one 
of the two wells of a double- well potential. What is then the rate with which the probability of finding the particle in 
this well decays in presence of an Ohmic-like environment? In addition, what are the asymptotic well populations? 
An additional manipulation of the potential barrier, i.e., a static bias or a time-dependent harmonic driving may be 
applied. In this work, we provide an analytic method to solve this non-trivial problem - also in presence of a time- 
dependent driving field - in a very general manner. We restrict ourselves neither to a large semi-classical potential 
barrier, nor to a weak system-bath interaction nor to weak driving fields. Our analysis is based on the real-time path 
integral technique which uses the Feynman- Vernon formulation as a starting basis. By treating the bath induced 
correlations between quantum paths within a generalized non-interacting cluster approximation, a generalized master 
equation for the diagonal elements of the reduced density matrix is derived. It turns out that the approximation is 
appropriate in the regime of moderate temperature and/or moderate system-bath coupling. A further simplification 
of the integro-differential equation leads to a Markovian approximated master equation whose rate coefficients are 
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obtained in the form of closed analytical expressions. By comparing the results of the full generalized master equation 
versus the Markovian approximated master equation and versus the numerical quasiadiabatic propagator path integral 
algorithm |21|, we conclude that the analytical approximation permits correct predictions for the decay process out 
of the initially populated potential well. The rate governing the long-time dynamics of the decay is obtained as the 
smallest eigenvalue of the matrix of the (time averaged) rate coefhcicnts. The dependence of this quantum relaxation 
rate and of the asymptotic population of the metastable well on the various physical parameters is investigated in 
detail. 

We stress that the developed method is not restricted to this specific problem but can be applied to many different 
other physical situations where a potential with a discrete energy spectrum can be assumed. A short summary of this 
present work has been published in Ref. |22|. 

Before we proceed, we motivate that the stated problem is not of formal academic nature but, in contrast, has 
several applications to real physical systems. For that purpose, we have collected numerous experimental works in 
the following subsection In the subsequent subsection [ B we briefly review the few existing theoretical works 
and discuss some of their shortcomings and inconsistencies which we attempt to overcome by our techniques. The 
rest of the paper is organized as follows: In Section |l| our specific model is introduced. The succeeding Section [II is 
devoted to the derivation of the dissipative real-time path integral which is cast in the discrete variable representation 
(DVR), i.e., the eigenbasis of that system operato r whic h couples to the reservoir. The example of the double-doublet 
system illustrates this transformation in Section III D , In Section [V we introduce an approximation to the so far 
exact real-time path integral expressions. This approximative treatment of the bath induced path correlations allows 
for the derivation of a generalized master equation (GME). This is shown in Section where also the lowest order 
expressions for the integral kernels of the GME are given. In Section VI we extract the leading rate for the decay out 
of one of the two potential wells. This is possible if one applies an additional Markovian approximation to the GME. 
A detailed study of the dependence of the quantum relaxation rate on the various model parameter is put forward 
Moreover, an investigation of the asym ptotic well population is presented. Finally, our conclusions 



in Section VIl 



together with an outlook are presented in Section VIII. 



A. Experiments 

Several experiments where dissipative multilevel systems are involved have been performed in many different physical 
systems. We report on four timely examples to motivate the importance and the need for a consistent and general 
theory for the above stated problem. 

The first set of experiments deals with quantum tunneling of magnetization in nanomagnets . A macroscopic 
sample of molecular magnets consists of a large number (typically lO"* — 10^^) of chemically identical magnetic clusters 
of the same magnetic size. They are regularly arranged on a crystal lattice. The single molecules have usually a large 
spin quantum number, typically S ~ 10. Experiments (see below) indicate a strong uniaxial magnetocrystalline 
anisotropy. It favors a doubly-degenerate spin alignment along the c-axis of the crystal, ms — ±5, and generates 
an energy barrier for the reversal of magnetization. This implies two-fold degenerate excited states corresponding to 
the spin-projections = ^{S — 1),±(5 — 2), ...,0 in a double-well potential |2j]. At sufficiently low temperatures, 
the spins can tunnel through the anisotropy barrier. Two such materials are currently studied in detail: The first is 
referred to as Mni2-acetate. It possesses a tunneling barrier of A[///cb — 62K (fee denotes the Boltzmann constant). 
Resonant tunneling of magnetization reveals itself as quantum steps in hysteresis loops which go along with maxima in 
the relaxation rate for specific values of an external magnetic field . The second candidate is known as Fcg and has 
the advantage that the anisotropy barrier is approximately three times smaller than in Mni2 (ACZ/fce — 22K). This 
property enhances the observed effect by several orders of magnitude as compared to the case with Mni2. For the Fes 
samples several experiments on quantum tunneling of magnetization have been reported as well p^ , p7{ . Especially 
interesting for us are the measurements by Wernsdorfer et al. ; those are performed at non-adiabatic driving fields 
and at temperatures where many doublets contribute to the dynamics. 



A second class of ex 
devices (SQUIDs) ||- 



3eriments addresses tunneling of the magnetic flux in superconducting quantum interference 
36| . The dynamics of the total flux threaded through the SQUID (or the phase difference 
across a current biased Josephson junction) obeys a collective motion of a macroscopic number of quasiparticles. The 
classical equation of motion for the flux dynamics maps to that of a particle moving dissipatively in a (symmetric) 
double-well potential. Its lowest left (right) well corresponds to one of the two fluxoid states (1) of the SQUID. For 
sufficiently low temperature, the transition between these states occurs via tunneling through the potential barrier. 
Measurements of the relaxation of a fluxoid state initially prepared in an rf-SQUID have addressed two different 
physical situations: The results in Ref. |^ have been interpreted as incoherent tunneling in a macroscopic two-state 
system and those in Ref. pOt] have been explained as resonant tunneling between two quasi-degenerate localized states 
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in different fiuxoid wells. The rate of tunneling out of the metastable well vs. the applied external flux exhibits a 
series of local maxima. These occur at those values of the external flux where the adiabatic energy levels of the biased 
SQUID potential form avoided level crossings. By applying a resonant time-dependent external rf-field, Han et al. ]3l| ] 
created a population inversion between the two adjacent fiuxoid wells. Furthermore, Silvestrini et al. reported 
the observation of energy level quantization in underdamped Josephson junctions above the crossover temperature 
which separates the classical from the quantum regime. Han et al. recently presented evidence for transitions 
between the fiuxoid wells due to cascaded, two-photon processes. In the latest work of this group, Friedman et al. 

report on the realization of a quantum superposition of macroscopic states in an rf-SQUID. Similar observations 
were recently made by the group of Mooij [ p5| , ^ where symmetric and antisymmetric quantum superpositions of 
macroscopic states of a dc-SQUID have been created. 

Another set of experiments concerns the tunneling dynamics of substitutional defects in solids [ p7{ . For instance, 
in a crystalline environment tunneling arises from defect ions which do not fit properly in the sites offered by the 
host lattice. The symmetry of the host crystal determines a complicated potential energy landscape with several 
degenerate minima for the defect ion. Golding et al. |3^] studied the relaxation rates of individual microscopic defects 
in a mesoscopic disordered Bi-metal. Since the sample dimensions were comparable to the phase-breaking length 
for quantum transport |^ at low temperatures, the sample's conductance was highly sensitive to the positions of 
the scattering centers. Their observations were found to be consistent with predictions of the dissipative two-level 
system at low temperatures. However, measurements at higher temperatures [ ^9| have indicated the failure of 
the two- level theory po[| . Furthermore, the study of thermally assisted tunneling of atomic hydrogen and deuterium 
in boron-doped crystalline silicon reveals that the relaxation rate calculated by a path integral centroid formalism 
differs from experimental measurements by two orders of magnitude. Finally, Enss and Hunklinger have pointed 
out several discrepancies between predictions of a semiclassical tunneling model, where the two wells are approximated 
by harmonic oscillators {harmonic-well approximation, see also Appendix and experimental measurements at low 
temperatures. They proposed an improved approach by taking into account elastic interactions among the tunneling 
systems to resolve these discrepancies. 

The last class of experiments concerns systems in chemical physics with the goal of controlling of chemical reactions 
pSj-^Sf. The hydrogen pair transfer in the hydrogen-bonded cyclic dimers of numerous carboxylic acids is used as 
a prototype system to study the relation between quantum tunneling and chemical kinetics. The measurements 
show that the free hydrogen-bonded dimers possess two energetically degenerate equilibrium configurations. They 
correspond to the two minima of a double-well potential. Both quantum tunneling and vibrational excitation are 
important for the transfer of the hydrogen pair. This has been studied experimentally in detail in Refs. |43j. A 
specific control scheme ("Hydrogen-Subway") has been proposed to steer intramolecular hydrogen transfer 

reactions in malonaldehyde by ultrashort laser pulses. The conventionally proposed method for the transfer consists 
in applying a laser pulse that lifts an initially localized wavepacket in the reactant region over the barrier thus allowing 
propagation towards the final product configuration. The new approach in Ref. p5[ | is to drive the wave packet not 
over but through the barrier. This is achieved by exciting higher lying doublets where tunneling occurs on a much 
shorter time-scale than in the lower doublets. The advantage of this new proposal is that it requires laser intensities 
which are considerably smaller than those used in the conventional approach. 

B. Prior theoretical approaches 

Previous theoretical works dealing with dissipative spatially continuous quantum systems, being driven or undriven, 
naturally fall into two classes: Approaches that are more of a numerical or analytical flavor, respectively. 

In Ref. 0, the harmonically driven double- well potential has been investigated numerically in presence of dissipation. 
For that purpose, a master equation for the reduced density matrix has been derived on the basis of the standard 
Born-Markov assumption [i^ . Subsequently, an analytical Floquet approach is used to derive the master equation. 
In doing so, an improved master equation has been obtained in Ref. | |47[ |. Here, the Floquet theory is applied on 
the level of the Schrodinger equation and the Born-Markov approximation is made for the quasienergy spectrum. 
In both cases, the system-bath coupling is treated perturbatively. This restricts the method to the weak- coupling 
regime. The same regime of a weak system-bath coupling was treated by Naundorf et al. p5[ . Also, standard Redfield 
(i.e., weak-coupling) techniques have been applied to derive a master equation. The specific shape of a laser pulse is 
determined in order to control hydrogen tunneling in a dissipative environment p5[ . In the strong coupling regime, 
the harmonically driven double-well potential has been studied in the context of quantum hysteresis and quantum 



stochastic resonance |11 . In this work, the system has been iterated numerically using the tensor multiplication 
scheme within the quasiadiabatic propagator path integral technique developed by Makri and Makarov [pT| . 

More analytical oriented works in the context of dissipative multilevel bistable systems have been performed by 
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several groups [^8|-p6| . 

The starting point in Refs. ||4^ , |49|] is a multilevel system with interdoublet transition terms (vibrational relaxation) 
which are not strictly derived from a continuous double-well potential; these are constructed phenomenologically. 
This leads to the assumption that the vibrational coupling occurs only between vibrational states located inside the 
same well. 

The group of Silbey |Q considered a static multilevel system. Additionally, only tunneling states differing by one 
quantum of vibrational excitation are assumed to be connected. Finally, it is assumed that the vibrational coupling 
within each well is the same for both wells. This a priori excludes the case with a static asymmetry of the potential. 

The group of Morillo and Cukier started out from a similar Hamiltonian like in [Q. They restricted the 
model further and included only the two doublets with the lowest energy, i.e., the so-called double- doublet system. 
The authors for the first time included a time-dependent driving which couples to a phenomenologically constructed 
dipole operator of the multilevel system. The system-bath interaction is treated perturbatively within a generalized 
Redfield approach. 

In a sequence of articles pO| , Dekker analyzed the real-time dynamics of a quantum particle in the dissipative static 
double- well potential ab initio by means of a multisite spin-hopping model. He derives the reduced quantum Liouville 
equation for the particle, thereby not restricting the dynamics to the lowest doublet only. The interdoublet vibrational 
dynamics is approximated by coarse-graining the density matrix elements on a time scale of many vibrational periods. 
It is further assumed that the localized states in the wells are approximated by the eigenfunctions of a harmonic 
oscillator {harmonic-well approximation). This latter assumption can be justified as long as the barrier height is 
large compared to the interdoublet energy gap. In this parameter regime, however, the application of the standard 
semiclassical rate theory [^j is appropriate, and even simpler to apply. In the deep quantum regime with low to 
intermediate barrier heights, this assumption increasingly becomes invalid and leads to considerable deviations of the 
approximated wave functions from the exact ones (cf. also Appendix ^). Also, the eigenenergies of the harmonic 
potential are considerably different from the exact ones for a shallow energy barrier. 

A related problem has been investigated in a series of theoretical works by Ovchinnikov and co-workers [^Tf^6| by 
applying semiclassical techniques. In Ref. [sijl Larkin and Ovchinnikov developed a method to calculate the decay 
rate of metastable voltage states of Josephson junctions. They constructed a kinetic equation for the probabilities of 
population of many energy levels. The transition probabilities are determined for a cubic potential in semiclassical 
approximation for weak system-bath coupling. This procedure assumes a decay into the continuum via quantum 
tunneling or thermal hopping. However, within confining potentials such as a double- well this assumption may be not 
justified. The effect of time-dependent driving is included within an approximation. The low temperature case where 
tunneling prevails is considered in Ref. for vortices moving in a washboard potential being weakly coupled to the 
environment. Also quasiclassical conditions have been assumed. The problem of divergent expressions for the decay 
rate at avoided level crossings is cured in Ref. fs^ ] where a two-level approximation at the avoided level crossings is 
invoked. The authors treat the problem within the harmonic well (i.e., quasiclassical) approximation for a constant 
spectral density of the bath modes, and for a weak system-bath coupling. The semiclassical expressions of Ref. ]5l| ] 
are applied to Josephson junctions (i) in Ref. ||54|] to calculate numerically the decay rate of the zero-voltage state 
for non-stationary conditions, and (ii) in Ref. | |55|] to study the influence of temperature for resonant macroscopic 
quantum tunneling. Finally, the theory is adapted to SQUIDs in Ref. to explain the experimental findings of Ref. 

. However, the theoretical results follow qualitatively those obtained from the standard WKB-approximation. The 
calculated decay rate differs from the experimental results by two-to-four orders of magnitude for small static potential 
asymmetries, i.e., with still large barriers, where the semiclassical treatment should yield rather good agreement. In 
contrast, for large bias asymmetries, one of the two barrier heights becomes rather small so that the semiclassical 
approximation is expected to yield worse results. The agreement with the experimental data turned out to be of the 
same order of magnitude. This inconsistency may be mainly traced back to the fact that the semiclassical treatment 
is not appropriate for a system in the deep quantum regime when only two to six levels lie below the energy barrier. 

In summary, no analytic treatment exists in the prior literature where tunneling and vibrational relaxation is 
investigated consistently in the regime where a finite number of discrete energy eigenstates rules the dissipative 
dynamics. This is so even for the situation that no time-dependent driving acts upon the system. While standard 
Redfield theory for a weak system-bath coupling is used frequently, the theory for the strong coupling regime for 
static as well as for driven multilevel systems is still in its infancy. The main objective of this work is to fill this gap 
in deriving analytical schemes that cover the physics in this prominent regime of a moderate-to-strong system-bath 
coupling. 
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II. THE DRIVEN DISSIPATIVE BISTABLE SYSTEM 



We consider a quantum particle with mass M, position operator q and momentum operator p moving in a one- 
dimensional double- well potential Vo(q) which may include a static asymmetry. The potential experiences a time- 
dependent external force, ssin(r2t), with field strength s and frequency Q. It is described by the Hamiltonian 

Hs(t) = Ho-q.ssinr!i= ^ + yo(q)-qssinm, (1) 

with 

being the asymmetric double-well potential. The quantity s denotes the static bias force. In absence of the asymmetry 
(e = 0), AU denotes the barrier height, and wq is the angular frequency of classical oscillations around the well minima. 

The energy spectrum of Hp follows from the time-independent Schrodinger equation with a static double-well 
potential Vb(q), i.e., Ho|n) — £„\n),n = 1, 2, ... . 

In absence of a static bias (e = 0) and for energies well below the barrier, the spectrum consists of a ladder of 
pairs of energy eigenstates (doublets) . The energy gaps within each doublet generally are several orders of magnitude 
smaller than the inter-doublet energy gaps and are responsible for the tunneling dynamics between the two wells. The 
large energy gaps are of the order of the harmonic oscillator energy gap Hloq associated with each well. For energies 
above the barrier, the energy gaps are also of the order of Hloo- Transitions between those largely separated energy 
eigenstates are termed vibrational relaxation. In presence of a static tilt (e 0), no general statement can be made. 
Spectra with typical avoided level crossings can occur as well as such with almost equally separated energy levels, cf. 
Fig. I a.). 

Following the common approach [p]-^|j5^] to model the influence of the environment by an ensemble of harmonic 
oscillators, the bath Hamiltonian Hb (including the interaction with the system) is given by 

j = l J -I J 

The whole system is thus described by the Hamiltonian H(<) = Hs(i) + Hb. In the case of a thermal equilibrium 
bath, its influence on the system is fully characterized by the spectral density 

jH^lx:^.(.-.,). (4) 

With the number M of harmonic oscillators approaching infinity, we arrive at a continuous spectral density. Through- 
out this work, we choose an Ohmic spectral density with an exponential cut-off, i.e., 

J{lo) = rjuj exp{—ijj juJc)- (5) 

Here, rj = with 7 being the strength of the coupling to the heat bath. Moreover, Wc ^ denotes a 

cut-off frequency being the largest frequency in the model. 

We choose a factorizing initial condition of Feynman- Vernon form [^7| . This means that at time t — to, the full 
density operator W(io) is given as a product of the initially prepared system density operator Pg(io) a-nd the canonical 
bath density operator at temperature T = l/k^P, i.e., 

W(to) - Psito) ^b' exp(-/?H^) , (6) 



r „2 



P ; 9 9 



where Zb — tr exp(— /SHg) and Hg = J2j=i 2 

In order to describe the dynamics of the system of interest we focus on the time evolution of the reduced density 
matrix. In position representation it reads 

p((7/,g};t) = tr,es(g/nja;,|U(t,io)W(io)U-i(i,to)|(z}H,4) , 

U(i,to) =Texp|-i/fi / Ii{t')dt'\ . (7) 



Here, T denotes the time ordering operator, W(fo) is the full density operator at the initial time to and tries indicates 
the partial trace over the harmonic bath oscillators Xj . 
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III. THE REDUCED DENSITY MATRIX IN THE DISCRETE VARIABLE REPRESENTATION DVR 



A. The Feynman- Vernon influence functional 

Due to our assumption of a factorizing initial condition in Eq. (|^) , the partial trace over the bath can be performed 
and the reduced density operator be recast according to Feynman and Vernon fs^ as 



P{<lf,q'f,t)^J dqo J dq'„ G{qf,q'f,t;qo,qQ,to)ps{qo,qa,to) , (8) 
with the propagator G given by 

G{qf,q'f,t;qo,q'o,to)^ Vq f^^ Vq' A[q]A*[q']TMq,q'] ■ (9) 

Here, A[q] — eKp{iSs[q]/fi} denotes the bare system amplitude, with Ss[q] being the classical action functional of 
the system variable q along a path q{t). J-'-pylq^q'] = exp(--(/)Fv['?, <z']/^) denotes the Feynman- Vernon influence 
functional. For our purpose, it is convenient to write the influence phase (t>Fv[q,q'] in terms of relative coordinates 
^{t') — q{t') — q'{t') and center of mass coordinates x(^') — li^') + respectively; it reads 



'/'Fv[x,c] ^ dt' dfimsit' - nan + imm' - nxin) 

J to J to 

+ / dt' {i{t' ) S it-t')+ licit' )Rit-t')} 

J to 

+ato){ms{t~to)~ [ dt'at')s{t' -to)} 



J to 

+^x(to){mRit-n- f dt'i{t')R{t' -to)}. (10) 

Jto 

Herein, S{t) and R{t) denote the real and imaginary part, respectively, of the bath correlation function Q(t), i.e., |^ 

Q(t) ^ S(t) + iR{t) = - [ duj^^^^ I coth^^il- COS ujt) + ismujt\ . (11) 
TTJo ^ I 2 J 

We evaluate in the following the reduced density matrix explicitly. It turns out that this is conveniently performed in 
the discrete eigenbasis of the position operator q. This representation is the so-termed discrete variable representation 
(DVR) The reason for this basis transformation is that only then can the influence phase, Eq. (p^), be evaluated 
at the eigenvalues q^ of q. This is shown in the subsequent section. 



B. Real-time paths in the DVR basis 

The time- independent double- well potential Vo(q), Eq. (||), possesses a discrete energy spectrum. The interesting 
temperature regime for us is that in which only a finite and small number of energy eigenstates is thermally significantly 
populated. A quantum mechanical description would not be necessary if the temperature is very large compared to the 
natural energy scale of the system. We assume furthermore that the time-dependent driving does not excite arbitrary 
high lying energy eigenstates of the static problem. Then, it is appropriate to consider only the M-dimensional 
Hilbert space spanned by the M lowest lying energy eigenstates of the static potential. The problem of a spatially 
continuous double- well potential is then reduced to a problem of a finite dimensional M-level system (MLS). The case 
of M — 2 (with e and s being sufficiently small) is the well-known (driven) spin-boson problem ||l|,||,D , while M = 4 
constitutes, for instance, the double-doublet system |^^. This reduction has been shown to be sensible for the case 
of the parametrically driven dissipative quantum harmonic oscillator js^ . There, the spatially continuous potential 
is appropriately described by a discrete M-level system with M = 3 to M = 6. 

Next we perform a basis transformation to the so-called discrete variable representation (DVR) Jsst . The new basis 
is chosen as the eigenbasis of that operator which couples the bare system to the harmonic bath. In our case this is 
the position operator q. We define the DVR basis {1?^)} according to 

{QnHq^) = qi^Sf^u , p,v=l,...,M. (12) 
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This basis follows from the energy eigenbasis {|m)} by inserting the identity I = X]m=i yielding \q^) = 

^^^j^(m|g^j)|m). This step allows to transform the description of the dynamics as transitions between energy eigen- 
states to a hopping among the M discrete position eigenvalues of the spatial grid. While for the static symmetric 
case, e = 0, the position eigenvalues are located symmetrically on the q-axis with respect to q = 0, this is no longer 
the case in presence of a static bias e ^ 0. 

To describe the dynamics in the DVR basis, we define a quantum mechanical path q{t') along which the system 
evolves in time. It starts out at time t' = to in the state q{t' = to) = 9^0 ^^"^ evolves via N jumps between the M 
discrete states into the final state q{t' = tff) — q^^. The full time interval is split into N short time intervals such that 
the jumps happen at times t' — tj. The intermediate states are labeled by q^^, where /Zj — 1, . . . , Af is the quantum 
state index, and j = 1, . . . , N — 1 denotes the time index. The full path is assumed to be a sequence of constant paths 
segments according to 

N-l 

q{t') = -q^^Mt' - ii) + E 1h - tj) - - i.+i)] + - %) ' (13) 

where Q{t) is the Heaviside function. Thus, upon switching to the center-of-mass and relative coordinates x(i') — 
q{t') + q'{t') and S,{t') = q{t') - q'{t'), respectively, (cf. Eq. (|l|) and Eq. (|l|) below), the double path integral over 
the M-state paths q{t') and q'{t') in Eq. is visualized as an integral over a single path that jumps between the 
Af ^ states of the reduced density matrix in the (q, g')-plane. The total number N of jumps is given by the sum of the 
number of jumps for the paths q and q' , i.e., N = N + N'. 

Fig. |l| illustrates this idea for a general Af-state system described by its M x Af density matrix. Two paths are 
depicted: one (full line) starts in the diagonal state {qi, q'l = qi) and jumps in TV' = 3 horizontal jumps and in = 2 
vertical jumps to the final diagonal state (93,93 = 93). It visits four intermediate off-diagonal states (filled circles). 
The second path (dashed line) starts in the diagonal state (92, 92 = ^2) and travels via two intermediate states to the 
final diagonal state (gA/, q'm ~ 1m)- 

The paths in the relative and center of mass coordinates read 

+ W-«0(i'-^Af)> (14) 

and 

X{t')^q{n + q'{t') 

N-l 

= -XMo.oe(t' - ti) + E ^A^.-. [®(*' - - ®(^' - ^j+i)] 

+ X,.«.«0(i'-ijv). (15) 
Herein, the path weights are given as 

^^•^^=1^-1'^, (16) 

and 

Xh'^j ^ 1h + I'l^j ■ (17) 

In this discrete notation, the index /i refers to the path q and the index v to the primed path q' . The time intervals 
in which the system is in a diagonal state of the reduced density matrix are called sojourns. They are characterized 
by ^{t') = and x{t') 7^ 0. The time spans in which the system is in an off-diagonal state are called clusters. The 
clusters are characterized by £,{t') 7^ and x{t') 7^ 0. This is different from the spin-boson problem [|l],|[|^ where the 
off-diagonal states (blips) are characterized by £^{t') 7^ and x(^') — 0- Upon determining the derivatives of the paths 
with respect to the time variable t', we find 

N 

c(0-E^^'5(t'-t,) (18) 
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and 



N 

x{t')^Y.^At'~h). (19) 

3 = 1 

Thereby, we have introduced new paths weights according to 

= ^tijV] ~ C/jj-ii^j-i (20) 

and 

(21) 

with j ~ 1,...,N. For j — 0, we define = Cmoi^o ^-^^ Xo = X^.o'^o- Hence, a path with N transitions at times 
ii, t2, . . . , t^r can be parametrized by two sets of path weights {xoj Xi, X2: ■ • • j Xw} and {CojCii'f2, ■ ■ • I'fjv}- In the 
influence functional the paths are coupled. The situation mimics the case of interacting electrical charges. Thus, the 
paths weights in Eqs. (po|), (|2l] ) are termed charges. In the discrete notation, the real-time path integral expression 
(^ assumes the form 

l-^{t)=^N rx{t)=XN 

= J2 2?xS[x,^]^Fv[x,C]Ppo.o- (22) 

Here, = ^^[qjy^* [g'], and the influence phase takes on the form 

N l-l N l-l 

0Fv[x, = - 51 E - - * E E - • (23) 

1=1 j=0 1 = 1 j=o 



C. The population of the left well 



Since we are interested in the decay of the population of one (metastable) well of the bistable potential, say the 
left well, we define the quantity of interest to be the sum of the populations of those L DVR-states jg^), fJ, = 1, L, 
which belong to the negative position eigenvalues q^j^, i.e., those which are located to the left from the zero. This yields 

L 

Picft(t) = E/'A"M(0- (24) 

In absence of a static bias, i.e., e = in Eq. (|l|), the energy eigenfunctions occur in pairs of symmetric and 
antisymmetric wave functions. This implies a choice for an even number M of states. Then, half of the position 
eigenvalues is on the left side and the other half is on the right side of the position point of reflection symmetry, being 
at g = for Vb(q) in Eq. (^. The consequence is that for the population P\cit{t) of the left well, usually L = M/2 
DVR-states are relevant. However, in the case of a finite static asymmetry, no such general statement can be made. 

To determine P\cit{t) in Eq. ( p^ we focus on the case that the final state (/iat, vn) of the system will be a diagonal 
state, i.e., 

VN = i^N ■ (25) 

Since then q{t) = q'{t), it follows that ^{t) = in Eq. (|l^. 

The initially localized wave packet is assumed to be a superposition of energy eigenstates . The transformation to 
the DVR-basis generates an initial system density matrix Pfj^g^a which generally is non-diagonal, i.e., 

1^0^ Ho . (26) 

Accordingly, we keep the general initial conditions Pfj,gua 

^ in Eq. (M). 
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We proceed to the expUcit evaluation of the path integral in Eq. (^2|) with the boundary conditions given in Eqs. 
(H) and 

To determine the transition amplitudes of the bare system we consider a discrete path starting in a general initial 
state that ends in a diagonal state. It is described by a sequence of pairs of state labels 

{^J.a,lya) (^2,2^2) ^ > i^J■N,l^N) = i^J'N,^J■N) ■ (27) 

The first symbol of each pair belongs to the horizontal direction and labels the rows of the reduced density matrix. 
The second symbol corresponds to the vertical direction and labels the columns. This implies that for a horizontal 
jump the first index remains constant, i.e., (/ij, Vj) — > (/ij+i, i^j+i) = (/ij, t^j+i), while for a vertical jump the second 
index is unchanged meaning {fjLj^Vj) (/ij+i, j/j+i) = 

We are interested in the probability amplitude of finding the system in state (/^j+i, Vj+i) after a time Ai = tj+i —tj 
having started from This quantity is given by the time evolution operator of the bare system. We find 

for a vertical jump, i.e., Vj+i = h'j, the amplitude (g^^-^j | exp{— iHoA</?i}|(7^^.) and for a horizontal jump, i.e., 
/ij+i = fij, {qt,.^-^\exp{+iH.oAt/h}\q^.), respectively. The relevant part of the system Hamiltonian Hs(i) in Eq. (|l|) 
is the time-independent part Hp since we are interested in the cases qp, ^i q^, and Qi^^^ 7^ Qu ■ Taking into account 
the exponential operator up to linear order in the argument, i.e., exp{±iHoAt/?i} w I± iHoAt/h, and using the 
orthogonality relation {qi\qrn) — ^im, the result for the transition probability amplitude per unit time At is obtained 
as ztiAj/2. Here, the factor of 1/2 is extracted to have the same convention as in the spin-Boson-problem. The 
factors Aj for a horizontal jump are defined according to 



and for a vertical jump 



A, =A,^.^,,^ =^(g,^.^JHo|rz.^), (28) 



A, = A^^.^,^^ = |(g^^^JHo|g^^) , (29) 



respectively. The + (— ) sign belongs to a horizontal (vertical) transition in the reduced density matrix. The different 
signs for horizontal and vertical direction reflect the fact that the bare transition amplitude A[q] belongs to the vertical 
direction, while the complex conjugate transition amplitude A* [q'] belongs to the horizontal direction of the reduced 
density matrix. 

The amplitude to stay in the j-th off-diagonal state lasting from tj to tj^i depends on the time-dependent 
diagonal elements of the bare system Hamiltonian in the DVR-basis. It is given by the so-called bias factor 

exp (i /^*'+^ dt'[E^^ {t') - E,^ (<')]) , where 

E^ii') = \{q^.^'i^s{t')W,) = ^(F^,. -9^,ssinm') (30) 

with F^. = {q^. |Ho|(7p^). For the entire evolution from to tot^, N oi these factors are multiplied, yielding the overall 
contribution exp{i YIi^^q It'*^ dt'[E^. {t') — E^. {t')]\- This defines the transition probability amplitudes of the bare 
system in a unique way. 

The functional integration over all continuous paths in Eq. (^) turns into a discrete sum over all possible path 
configurations {fJ-ji^j} in the DVR basis and an integration over all intermediate times {tj}. In formal terms this 
implies 



jvijvxr^j V{t,} , (31) 

where we have introduced a compact notation according to 



V{tj}= I dtN / dtN^i--- I dt2 / dti (32) 

tc 'J to ^to J to 'J to 



for the time ordered integration over the N transition times tj in Eq. (22). 

Collecting all parts we obtain the dissipative real-time path integral for the diagonal elements of the reduced density 
matrix of an M- level system in the DVR-basis, i.e., 
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M oo „t f JV-1 „tj+i 

„ , ,„ — 1 7V7 1 •'in f. .. 1 I A — n ^ t-i 



Af-1 / , \ N 

X 

{AT i-1 TV /-I 1 

E E ^'^(*' - + ^ E E - f • (33) 
/=1 j=Q 1=1 j=0 J 

In this expression, the sum over aU possible path configurations {nji^j} in the spirit of Eq. ( p7|) has to be performed 
with Sj = 0(1) for a horizontal (vertical) jump. 

Several comments on this quite comprehensive path integral expression are apposite: First, the path integral in Eq. 
( p3| ) is given in its most general form and is formally exact because no approximations, neither on the form of the 
system Hamiltonian nor on the type of the system-bath interaction, are made. This method could be applied to any 
problem where a potential with a discrete spectrum is given, and where the coupling to the heat bath is mediated via 
the position operator. The main ingredients are the matrix elements of the system Hamiltonian, being represented 
in the DVR-basis, and the position eigenvalues via the paths weights. No specific requirements on the shape of the 
external driving have been made; even a stochastic driving force (such as multiplicative noise) can be included. 

In the case of only two levels, i.e., M = 2, Eq. (^3|) reduces to the well-known expression for the (driven) spin-boson 
problem ||l|,§,||. There, the problem simplifies due to the fact that the path weights during the time evolution take on 
only two values, corresponding to the two states localized in the left and in the right well of the potential. This means 
that the path flips between a sojourn and a blip at each jump. This implies that the spin-boson path integral assumes 
the form of a power series in the tunneling splitting Af = £2 — Si of the two lowest levels. This is not necessarily the 
case for a general M-level system where a path can travel around, visiting many off-diagonal states, before ending in 
a diagonal state. Certainly, such a path becomes less likely the longer it remains off-diagonal. This is due to damping. 

The path integral is not tractable in its most general form without assuming further approximations. Such an 



approximation is developed in the following Section IV. However, to gain insight into the physics behind the formal 
expression (^3|), we introduce in Section HID the example of the so-termed double- doublet system and discuss the 
transformation to the DVR-basis. It refers to the case where two doublets in a symmetric double-well potential, Eq. 
(0), are localized below the barrier, i.e., the case M — 4. 



D. An example: The symmetric double-doublet system 



We illustrate the general method with the example of two doublets below the barrier in the double-well potential, 
Eq. (^. Choosing Af = 4 generates the first non-trivial extension to the familiar spin-boson problem. 

For the sake of simplicity, but without loss of generality, we consider the symmetric potential, i.e., we set e = 
in Eq. (|l|). For the isolated system the energy spectrum follows from the time- independent Schrodinger equation as 
Ho|n) = £„|ri), n = 1, 2, .... The two lowest doublets ?iAf = £2 — £i, and fi,Af = £4 — £3 are separated by the energy 
gap fuIjQ = ^{£4 + £3) — ^{£2 + £1) ^ . The interdoublet frequency ljq is of the order of the classical oscillation 
frequency uq and becomes equal to it in the limit of high barriers when the two intrawell oscillators approach harmonic 
oscillator potentials. With the objective of the decay of a localized state in mind, we start from the so-called localized 
basis. It is this basis which is favorably used to describe the tunneling dynamics. It follows from the energy eigenbasis 
by a unitary transformation according to 

\Li) = -i=(|l)-|2)), |i?i> = -i=(|l) + |2)), (34) 
\L2) = ^(|3)-|4)), \R2) = -i=(|3) + |4)). 

These states are localized in the left {\Lj)) and in the right (\Rj)) well with lower (j — 1) and higher {j — 2) energy, 
respectively. The localized states are depicted in Fig. ^ a.) in position space. Shown is the double- well potential 
(thick solid line) for a barrier height of Eb = AUJhuJo = 1.4 (we use in the figures dimensionless quantities according 
to the standard scaling defined in the Appendix \A\ Eq. ( [Al| ) ). The energy eigenvalues £1, ...,£4 are marked by thin 
solid horizontal lines. The wave functions (gjii) (solid line) and {q\L2) (dashed-dotted line) are localized in the left 
well, and the wave functions {q\Ri) (dashed line) and {q\R2) (long dashed line) are localized in the right well. In the 
literature |Q , these localized states in Eq. (^) are sometimes approximated by the eigenstates of harmonic potentials 
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shifted to the position of the well minima, cf. Appendix This approximation is justified for large barrier heights 
where, however, semiclassical techniques l^] to determine the quantum relaxation rate are already applicable. By use 
of basic algebra, the matrix for the bare system Hamiltonian of the double-doublet system in the localized basis is 
calculated to be 

Hl5&S= E -^m{L^\ + \U){R^\) 
1=1,2 

+hnJo{\R2){R2\ + \L2){L2\), (35) 

with frequencies Af and uJo defined above. The position operator in this localized representation then reads 

q'°'= E a.,(|i?.)(i?,|-|L,)(L,|) + 6(|Li)(i?2| 

+ \R2){Li\-\Ri){L2\-\L2){Ri\), (36) 

where an = (l|q|2), 022 = (3|q|4), ai2 = flai = ((l|q|4) + (2|q|3))/2 and b = ((l|q|4) - (2|q|3>)/2 < Oy. Note that, 
in clear contrast to the spin-boson case M — 2, the position operator in the localized basis is nondiagonal. Since the 
energies in the Hamiltonian are of different orders of magnitude, i.e., hAf <C fi-Af <C tiOJo, the general time evolution 
of an initial state proceeds on different time scales. The coherent dynamics exhibits transitions between the wells due 
to tunneling. It occurs in the lower doublet on a time scale (Af and in the upper doublet on a much shorter time 
scale (Af )~^, being still long compared to the time scale uJ^^ of the interdoublet dynamics. The coupling to the heat 
bath is mediated by the position operator while the interdoublet transitions are responsible for vibrational relaxation. 



For the following analytical treatment, we simplify the approach by setting 6 = in Eq. (36). This is for the sake of 
an illustrative purpose only and has no impact on the path integral formalism introduced above. For specific results, 
the diagonalization of the position operator is performed numerically on the computer with 6 7^ 0. By means of 



ordinary diagonalization performed for the matrix in Eq. (36) the DVR-states read 



|ai) = v{\Li) ~ u\L2)) , \l3i) = v{\Ri) - u\R2)) , (37) 
1^2) = v{u\Li) + \L2)) , m = v{u\Ri) + \R2)) , 



with \aj) being localized in the left (right) well, respectively. Here, v — l/vT+t? and u = (an + qa-i)/ai 

— {0,22 + 9Q2)/ai2, and (/q.. = — g/j^ denote the position eigenvalues: 



-(an + a22) T \ (an - 022)^ + 4a' 



/2. (38) 



The four DVR-states are depicted in Fig. || b.) for a barrier height of Eb — AU/huJo — 1-4, i.e., {q\ai) (solid line), 
{q\a2) (dashed line), {q\(32) (dashed-dotted Hue), and {q\(3i) (long-dashed line). On the g-axis, the exact eigenvalues 
Ou are marked by crosses (the eigenvalues are obtained by numerical diagonalization of the position operator in Eq. 
(|36|)). As expected, the DVR-states are localized around their corresponding position eigenvalue (7^. 

It is suggestive to call transitions between the left and right well, i.e., between \ai) and \f3j) as DVR-tunneling . 
These are characterized by the effective tunneling matrix elements 

Aa^.fJ, = i;2(Af + U^Af ) , Aa,,f3, = v\u^Af + Af ) , 

A^,p, = A^,p, = v^u{Af - Af ) , (39) 

which constitute a linear combination of the bare tunneling splittings Af and Af . On the other hand, transitions 
within one well, i.e., between \ai) and \aj) and between \f3i) and \Pj), may be termed DVR- vibrational relaxation. 
Those can be characterized by the transition matrix elements 

Aa^a2 = ^a2ai = ^Pif32 = ^^2^1 = = 2u^uZUo . (40) 

Due to parity symmetry, they assume equal values. The Hamiltonian of the double- doublet system in the DVR-basis 
can thus be written as 

HEi5t = - E ^?iA„.^^.(K)(/3,| + |/?,)(a.|)-i;iARR 

+ E ft(i^aJa.)(a.|+Fft|/?,)(A|), (41) 
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with Fa-^ = Fp-^ = v^v^TDq, Fa^ = = v'^ujQ. The operator R accounts for DVR-vibrational relaxation, i.e., 
R = |ai)(a2| + |Q!2)(ai| + \Pi){(32\ + \P2){Pi\- Thereby, the time-independent problem is fully characterized. The 
time-dependent driving s{t) = s sm{^lt) couples to the position operator q. The total system Hamiltonian in the 
DVR-basis thus reads 

n^^\t) = HE^f - s sin{nt) g„.(|a,)(a,| - • (42) 

1=1,2 

Note that in the DVR-basis the time-dependence enters only in the diagonal elements of the Hamiltoni an. So far, we 
have provided all required parameters for the path integral formalism developed in the previous Section [II C| , namely 



the matrix elements in Eq. ( p| ) of the Hamiltonian in the DVR-basis, and the position operator eigenvalues in Eq. 
(^). In the following Section, we return to the dissipative real-time path integral formalism and develop a suitable 
approximation scheme to the exact expression in Eq. (p3|). 

IV. THE GENERALIZED NON-INTERACTING CLUSTER APPROXIMATION 

In the context of dissipative real-time path integrals, a common strategy of approximative treatment is as follows: 
It concerns the treatment of the interactions between different paths which are induced by the coupling to the heat 
bath and which are described by the influence functional. The working idea behind the strategy is to neglect 
some of those correlations in order to get tractable expressions. One possible approximation within the spin-boson 
problem {M — 2) is the so-termed non-interacting blip approximation (NIB A) There, the interactions between 
off-diagonal states (blips) are neglected. In the NIBA the sojourn-blip interactions are disregarded except neighboring 
ones, and even those are treated approximately. Within the NIBA, the influence phase simplifies drastically and the 
path integral series in Eq. ( ^3| ) reduces to terms which are of lowest order in the level splitting Af . The NIBA 
can be justified in the case of Ohmic damping for high enough temperatures and/or large dissipation strength. In 
this regime, the average blip length is small compared to the average sojourn length, and the blip-blip and the blip- 
sojourn interactions can consequently be neglected. In fact, for finite temperatures and Ohmic damping, long blips 
are exponentially suppressed by the intrablip interactions. The NIBA fails for systems with asymmetry when the 
friction becomes weak; however, it becomes a systematic weak-damping approximation down to zero temperature for 
the case of a symmetric spin-boson system ■ 

Improved approximations take into account some of the correlations between the blips. One such step of an improved 
approximation has been denoted as the interacting blip chain approximation (IBCA) Q. There, the interactions 
of all nearest-neighbor blip pairs and the full interactions of the nearest-neighbor sojourn-blip pairs are taken into 
account in addition. This improved approximation confirms the validity of the NIBA in the stated parameter regime. 
It is valid also in an extended parameter regime where the NIBA already breaks down. 

The NIBA being applicable for a spin-boson system, i.e., a system with two tight binding sites has been generalized 
to the case with arbitrary many tight binding sites by Egger, Mak and Weiss ]6l| ]. In their work only tunneling 
transitions between nearest-neighbor sites are considered. The multisite paths along the discrete states of the reduced 
density matrix result in a sequence of sojourns (time intervals with the system being in a diagonal state) and clusters 
(time intervals with the system being in an off-diagonal state). It turns out that the corresponding path weights 
of the clusters sum up to zero. Consequently, the clusters can be considered as neutral objects. This suggests to 
neglect all interactions between the clusters yielding the non-interacting cluster approximation (NICA). In the time- 
mdependent problem considered in Ref. [Q, the time-integrations over the sojourn times in the path integral appear 
as convolutions. This feature makes the expression solvable by means of Laplace transforms. 

Motivated by the NICA, we here generalize it to the case with many levels by observing that the multilevel 
problem can be mapped onto the multisite one when also tunneling between non-nearest neighbor sites are considered. 
Moreover, we generalize the approach in Ref. [ |6l| by taking into account a time-dependent system Hamiltonian, 
together with a general initial reduced density matrix. 

The key argument in Ref. ||6^ refers to the overall neutrality of a cluster because the cumulative charge is zero. 
This is also the case for a general multi-level path integral in Eq. (|3^). To show this, we consider a general cluster at 
a diagonal state {fJ,k,Vk = Mfc) at time tk- It subsequently travels around among arbitrary many off-diagonal states 
and re-enters at time ti a diagonal state (fii , t/; = l^-i)- The cumulative charge Wd of this cluster is the sum over all 
individual path weights defined in Eq. (|2C|), i.e.. 
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j=k+i 

= , (43) 

because — [ih and vi = [ii. In this language, any path is just a sequence of sojourns and chistcrs, where the 
^-charges within each cluster sum up to zero. 

In general, the influence functional in Eq. ( p3| ) couples the x-charges inside each cluster with the ^-charges inside 
all other clusters. Similarly, all ^-charges are coupled to each other. These interactions consequently render the path 
summation intractable. 

Since the entire cluster can be seen as a neutral object which is only weakly interacting with all other clusters 
which are themselves neutral, it is suggestive to neglect all the intercluster interactions in the influence phase in Eq. 
(^). However, all, the mtracluster interactions, as well as all interactions of a cluster with the preceding sojourn are 
fully taken into account (see below). For a path starting in an off-diagonal state, we call that part of a path which 
precedes the first sojourn a semi- cluster. Within our approximative description, we neglect the interaction of this 
semi-cluster with all the later clusters but take into account the correlation of the first sojourn with the preceding 
semi-cluster. This "coarse-graining" is performed for general transitions. We call this the generalized non-interacting 
cluster approximation (gNICA). 

Before we exploit the consequences of this approximation, we discuss its regime of validity. The gNICA is justified 
when the average cluster length is small compared to the average sojourn length. This is fulfilled for high temperatures 
and/or strong damping. Far excursions from the diagonal state are damped exponentially, see Eq. ( [l0| ) for the influence 
phase. As such, the gNICA becomes exact in the limit of high temperatures. 

For the case of the spin-boson-system at low temperature T, small friction 7 and no bias (e = 0), the interblip 
correlations are only of second order in the coupling strength 7 while the mtrablip ones are of linear order in 7. Hence, 
the gNICA is a good approximation down to zero temperature |^. However, with M > 2, lowest-order contributions 
to interblip correlations arise due to the non-zero diagonal elements F^, in Eq. (|l|). This yields a rough condition 
of validity for the gNICA; it reads: HA^^^ < k^T or A^^^^ < 7, where A^^^^ = max{Af , Af , . . . }. On the other 
hand, an upper limit for the allowed values of the damping constant can be extracted by the following argumentation: 
The damping leads to a level broadening of the unperturbed eigenenergies. This is seen best in the form in Eq. ( p^ ) 
for the influence phase. The damping could be viewed as an additional contribution to the bare system propagator. 
The contribution is of stochastic nature and implies the level broadening. In order that a tunneling description 
makes sense, this frictional level broadening should not exceed the bare interdoublet level spacing, i.e., 7 ^ t^o- This 
condition is not really restrictive because friction strengths of the order of the oscillator frequency, i.e., 7 « wq, would 
indeed strongly suppress quantum effects. 

V. THE GENERALIZED MASTER EQUATION IN THE DISCRETE VARIABLE REPRESENTATION 

A. General derivation 

First, we address the non-driven case. We start by observing that every path which begins and ends in a diagonal 
state can be seen as a sequence of p clusters punctuated by sojourns. For paths starting out at an off-diagonal 
and ending in a diagonal state, also the initial semi-cluster appears. Within the gNICA prescription, it is now 
straightforward to see that the integrations over the sojourn times in Eq. ( p3| ) appear as convolutions! To use this 
property effectively, we switch to the Laplace transform Pp,kp.nW- It then follows that the integration over each 
sojourn contributes a factor A~^, while each cluster yields a factor which depends on the number of charges and on 
their configuration inside that particular cluster according to Eq. (^3|) . This very point is elucidated with an example 
in Appendix there, we present in detail the contribution of one specific path to the full path sum. In a second 
step, we generalize this idea. 

We consider transitions from the initial state (/ig, I'o) at time ta to the final state {(J-n, (J-n) at time t^- In doing so, 
we must distinguish between two cases: (i) the initial state is a diagonal state, i.e., hq = vq and (ii) the initial state 
is an off-diagonal state, i.e., /^o 7^ t'o- We separate the contributions to the path sum and obtain then for the Laplace 
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transform A) = /j, dte ''^Vmiva'n (0 the expression 

M M 



(A), (44) 



where p^„^jv,d (P/JN/ijv.o) denotes the contribution of the diagonal (off-diagonal) initial part, respectively. 

To proceed, we need to consider an arbitrary cluster which begins in the diagonal state (/Zj, Vi = fii) at time ti and 
ends in the diagonal state (/ij, Vj = ^ij) at time tj. We sum over all thepath configurations and denote this collected 
contribution the cluster function h^.^.{X). Proceeding as in Appendix |D| yields for the cluster function 



roc POO 

/iMjM.(A) = E / d,Ti- ■ ■ c?Tm_i exp{-A(ri H \- t„i-i)} 

m=|i-j| 

Eexp dt'[E,,^Jt')-E^,^,{t')]\ 



{^fci-fc} k k=l 



fc=0 ^ ^ 

E E ^'+*'^ ( ^ ) + * E E ^'+»-^ ( ^ ) 

1=1 k=l \n=k ) 1=2 k=l \n=k ) 



(45) 



with the difference times = tfe — ifc-i and with the conventions and notations taken from Eq. (paj). 

Each contribution to /3;^„;i„.D(A) can be viewed as a sequence of sojourns punctuated by clusters. Thus, in the first 
case (i), we sum up the contributions of all paths which start in (/^oj^'o = Po) ^^nd end in {untI^n = Pn) and which 
contain p clusters starting in some intermediate diagonal states (crfc,(Tf.) and ending in (cTfe+i, Cfc+i), i.e., 

pSm„,d(A) = E ^/»-i.Po(A)^/»..,..(A) . . . /i^„,.,(A)i , (46) 

where the sum runs over all possible intermediate diagonal states a — 1, . . . , M . The factors 1/A are the results of 
the integration over the sojourns, see Appendix 

In the second case (ii), where the initial state is an off-diagonal state, we assume that the path travels among 
off-diagonal ones and hits after d transitions for the first time a diagonal state [kj,, Kd) at time td- This part of the 
path is termed semi-cluster and the interaction with all the other clusters is neglected according to the gNICA. The 
sum of all such semi-clusters that start in (/io, t'o) and end in (k^, k^) results in a semi-cluster function fK^K^.tioi^oW- 
From time td on, the formalism from (i) is applied. Summing over all possible diagonal states Kd = 1, . . . ,M yields 
the contribution to the path sum 



,,o(A) = E { ■l'i^dKd,fJ.oi^aW E ^'*o-i.Kd(A)^'l(T2,(Ti('^) • ■ • ^A'iv,o-p(-^)-^ / 

Kd = l CTi,Cr2,...,(Tp J 



P?iM«,o(A) = E { f^d^d,,oM E 4'^-i.«d(A)T/^.„., (A)... (A)4 \ , (47) 
where the semi-cluster function is given as 
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+ --- + r„_i)} 



J2 '^xp V II 



dt'[E^,{t')~E,,{t')\ 



y k=i ■^Si=i '^i 

fe=i ^ 

( m l-l 

EE^'^(E-"U^+^EE^'^(E-" 



'i-i 



'/-I 



Xfc 



\n—k 



1=2 k=l 



\n=k 



(48) 



Defining the cluster matrix 7i(A) witli tlie matrix elements /i^^.^. (A), we can rewrite Eq. (46) and the inner sum in 
Eq. ( ^ ) as a matrix product, i.e., 



n{\) 



(p) / 1 



7i(A) 



A 



(49) 
(50) 



In a last step, the summation over all possible numbers p of clusters within a path has to be performed. This last 
sum can be formally recast, yielding 



(A) 



1 



A-7^(A) 



PfJ-NfJ-N -,0 



r 1 

(A) — y ] I /KtjKd.^ioi'o (A) j \ _n- 



-H{\) 



(51) 
(52) 



We insert this result into Eq. ([44|), exchange the order of summation in the second term of the r.h.s. and end up with 

M f 1 1 *^ 



(A) - E ^A^opo I ^ _ I ^ I A - H(A) } 



«Kd«:d(A) 



(53) 



with 



A/ 



(A)= E 

(A), 



(54) 



Eq. ( |53| ) can be viewed as a vector equation with two vector-matrix products on the r.h.s.. For convenience we 
introduce a vector-matrix notation. (A) then appear as elements of a vector p{\). The initial populations P^q^q 

arc arranged in the vector po and the initial off-diagonal elements are contained in the vector /(A) with the elements 
*K£iKd(A). In this notation, Eq. ( ^3|) reads 



P(A) 



A-7i(A)^° ' A-7^(A)^ 
Multiplying Eq. (|5^) with the inverse of the matrix and rearranging the equation, we find 

A/?(A) - po = W(A)p(A) + /(A) . 
Finally, we perform the inverse Laplace transform and end up with the equation 

'p{t) = f dt' Hit - t')p{t') + I{t - to) 

J to 



(55) 



(56) 



(57) 
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or in the original notation for the single components 



M t 

= / dt'H^,{t-t')p,,{t') + I^{t--to), M=1,...,M. (58) 

The overdot denotes the derivative with respect to time t. The initial conditions for Eq. ( |58| ) are p^,y{to) — p^gi^Q- 
Eqs. ( |57| ) and ( p8| ) are of convolutive form since a time-mdependent Hamiltonian was assumed. A similar, although 
technically more involved line of reasoning must be used for the driven case. We find equations similar to Eqs. ( |57| ) 
and ( p8| ) with Hf_i^{t — t') r\ Tifii,{t,t') and /^(t — Iq) rv 7^(t,fo). To be explicit, the elements of the rate matrix 
T^iiuit, t') are in the general time-dependent case given as 



^1 + ! 

dt"[ii;^^(t")-i?.,(t")] 



Af-l / • \ W 

X 



j=0 



N l~l N l-l 



X 



i E E - + * E E ^'^(^' - ) ■ (59) 

[ ;=i j=o 1=1 j=a 



The inhomogeneity /^(i, io) arises because of the contributions of the non-diagonal initial states; its explicit form 
reads 

M 

= E ^A^o-oE / ^{^^} E e^PV E / rfi"[i?M.(i")-i?..(i" 



m— 1 / ■ \ m 



i-l m /-I 



«^p i E E - + * E E ^'^(*' - r • 

[ 1=1 j=0 1=1 j=0 



The integro-differential equation ( p8[ ) is called the generalized master equation ( GME); it constitutes one central result 
of this work. 

In the following, we will see that the inhomogeneity in Eq. ( |60| ) plays an important role at short times. However, 
it will become exponentially suppressed at long times reflecting the fact that the asymptotic state is independent of 
the initial preparation. 

We note that this integro-differential equation ( |58|) is represented in the DVR-basis for the diagonal elements of 
the reduced density matrix p(t). For all practical calculations, the kernels in Eq. ( p9[ ) and the inhomogeneities in Eq. 
( |60| ) have to be determined up to a certain order ©(A^). In practice, this means that = oo as the upper limit of 
the summations has to be replaced by a finite value. 

Some comments to elucidate the physical content of the GME (^8|) are in order: The transformation of the problem 
from the localized basis to the DVR-basis maps the dynamics of the particle in the spatially continuous potential onto 
a hopping process of the particle on a spatially discrete grid. The grid points are the discrete positions characterized 
by the eigenvalues g„, k = 1, Af, of the position operator q, according to Eq. ( p^ ) 



Next, we consider the example of the double-doublet system which has already been introduced in Section [II D. 



We give the explicit expressions for the kernels in Eq. (p3) in the GME up to second order in Aj and illustrate the 
damping mechanism further. 
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B. The leading order approximation 



In this Section we investigate the GME, Eq. (|5^), with the kernels in Eq. ( |59| ) and the inhomogeneities in Eq. ( |60| ) 
derived to lowest order in the Hamiltonian matrix elements Aj. To illustrate the general scheme, we describe the 
method for the case of the double-doublet system with AI — 4 levels. 

As we want to evaluate the relaxation rate of an initially localized wave packet in one of the wells, say, the left well, 
we prepare the system in an equally weighted superposition of symmetric and antisymmetric wave function belonging 
to the lowest doublet, i.e., 

p{to) = \L,){L,\, (61) 

where \Li) = -^(|1) — |2)) and \n),n — 1,2 are the nearly degenerate energy eigenstates of the static, symmetric 
Hamiltonian Hp (cf. Eq. (||) and Eqs. (p^)). Transforming this initial state to the DVR-basis via Eq. (^^, we find 

p{to) = w^(|ai)(ai| + M^|ck2)(a2| + u\ai){a2\ + -u|a2)(ai|) (62) 

with the parameters u and v defined below Eq. ( ^7| ) . We note that the initially prepared localized state is characterized 
by a nondiagonal density matrix in the DVR basis. The diagonal elements in Eq. (^2|) enter as initial conditions for 
the first part of the r.h.s. of Eq. (|5^), while the off-diagonal elements determine the inhomogeneity. 

To first order, i.e., with one jump, no transition from an initial diagonal state to a final diagonal state is possible. 
To achieve this, at least two jumps are necessary. However, transitions starting in an off-diagonal state and ending 
in a diagonal state are possible within one jump. This means that a first order contribution appears only in the 
inhomogeneity of Eq. (|58|). The relevant transitions are the jumps ending in the diagonal state {ai,ai), i.e., {ai,a2) 
(ai, ai) and (a2, ai) ^ Q;i); and the jumps ending in (a2, 02), i.e., (ai, 0:2) (a2, 0:2) and (02, ai) — > (a2, 0:2), 
respectively. From Eq. (BSl) , it follows that each path traveling "above" the diagonal has a corresponding mirror path 
traveling "below" the diagonal. The mirror path yields a contribution to the path sum being the complex conjugate 
of the upper path. Using this feature and the fact that PqiqjI^o) = Pa2aii^o) = uv'^, we obtain for the inhomogeneity 
the following first-order expression 

lli'Kt^ to) = (^^iQi - 3^02) uv^ AqjQ2 exp {-(qttj - (JqJ^ S{t - to)} 

X sinj^ dt' [S„,(i') - ^ai(t')] - (fci - fcj' R{t - io)| . (63) 

Here, A^, = (/i|HQp^|j/), ^jl ^ v, are the off-diagonal matrix elements of the system Hamiltonian in Eg. (p]), see also 
Eq. (flOj). Note that only vibrational transitions within the initially populated well contribute in Eq. (|63[). Moreover, 
the (7„ (k = Q;i,a2) denote the position eigenvalues, see Eq. (^8|), and E^(t') = — QkS sin(Sli') are the time- 
dependent diagonal elements of the system Hamiltonian, see in above Eqs. (^), (^) and (^). The influence of the 
bath enters via the real and imaginary part of the twice integrated bath correlation function, i.e., S{t) and R(t), 
respectively, see Eq. (|ll|) and Appendix The conservation of probability is reflected with the opposite signs of the 
Kronecker symbols d,^K- From Eq. (^3|) it clearly follows that the contribution of the initial off-diagonal states are 
damped exponentially on a time scale determined by the damping constant 7 and the temperature T. We recall that 
the lowest order of the contribution of the integral part of the GME, Eq. (p8|), is of second order. This implies that 
the contribution of second order to the inhomogeneity should also be taken into account for a consistent treatment. 
However, we refrain from writing down the complicated second order term which would yield only minor physical 
insight. It can be neglected anyhow when investigating the long-time dynamics in the following Sections. 

The lowest order for the kernels in the integral part of the GME (^sj) is the second order, because at least two jumps 
are required starting in a diagonal state to end again in a diagonal state. We use once more the feature that each 
path traveling above the diagonal has a mirror path traveling below the diagonal, yielding the complex conjugate of 
the upper path contribution. We then obtain for the GME kernels the leading order results 

W(fj(t, t') = ^ exp {-{q^ - q,)' S{t ~ t')} 

X cos dt" \E,{t!') - S^(t")] - (g^ - qufm - t')| , tii^v. 

(64) 
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The conservation of probability implies for the diagonal kernels the condition 



M 



n^^J{t,t')^~yn';^J{tX}- (65) 



We emphasize here that the lowest-order expression in Eq. ( |64| ) is applicable to a general number M of levels. The 
explicit example of the double-doublet system with Af = 4 is used for illustrative purpose only. 

Note that the Tijjiu represents the transition probability for a path starting in the diagonal state {v, v), then jumping 
to the off-diagonal state {v, fj,)/{fi, v), and finally ending in the diagonal state (p,^). In clear contrast to Eq. (|6^), 
now tunneling and vibrational relaxation both contribute in Eq. (|64|). 

The structure of the GME with the kernels H-iiy restricted to leading, i.e., second, order is similar to that one 
obtained for the driven spin-boson system within the non- interacting blip approximation (NIB A) and to that 

one for the dissipative tight-binding model within the non-interacting cluster approximation (NICA) performed to 
lowest order The main difference to these GMEs is that in our case the factors {q^ — q^)"^ enter as prefactors 
for the damping constant 7 in S(t) and R{t), respectively. Since they arise from non-nearest neighbor hopping on 
a non-equally spaced grid of DVR eigenvalues, they are not equal for all transitions. This means that transitions 
between far away lying DVR-states are stronger damped and therefore less probable compared to those lying close to 
each other. This insight is especially relevant for the tunneling transitions from one well to the adjacent. Then, the 
main contribution to the dynamics comes from those two DVR-states which lie closest to the barrier within each well. 

One remark on the notation should be made: In the following, we use the superscript in the kernels Tifii, when 
they are utilized in second order. Whenever this superscript is omitted, the respective formula is valid to any order 

of H^i.. 

The generalized master equation (^) is clearly not solvable analytically in closed form, not even with the kernels and 
the inhomogeneities approximated to lowest order. Thus, in Appendix O, we provide a numerical iteration algorithm 
to obtain a numerical solution. 



C. Comparison with numerical ab-mitio path integral simulations 

In this subsection we compare the results for Pioft(i) obtained from the numerical solution of the GME js^ ) with 
those of the numerical iterative algorithm using the method of the quasiadiabatic propagator path integral QUAPI 
of Makri ||2l| . It is known that the QUAPI technique yields reliable results for time-dependent spatially continuous 
confining potentials |5^. Hence, we use it here as a reference in order to check the gNICA. 

We present results for the double-doublet system M = 4. Fig. |^ depicts the outcome for Picft(i) for the symmetric 
{e = 0) and for the asymmetric (e = 0.08) system. Each figure contains three lines: (i) the results of the full generalized 
master equation (full line), (ii) findings of the QUAPI algorithm which are used as a reference, and (iii) the outcome 
of a Markovian master equation which is introduced in the following Section We postpone the entire discussion 
of the Markovian results to the following Section. We find a very good agreement, both for the symmetric as well as 
the asymmetric system. We note that for the asymmetric case, the full GME is solved onlyup to t — 1000 due to 
the necessary choice of a very small At = 5 x 10~^ (for a detailed discussion see Appendix |c|) . The QUAPI results 
have been obtained with K = A (the number of memory time steps, see Refs. pl| , p9{ for details) and At = 0.1 for the 
symmetric and At = 0.35 for the asymmetric case. 

The same very good agreement is found for the case with resonant driving, i.e., s = 1.0,57 = ZJq = 0.815 which 
is depicted in Fig. ^ a.) for T 0.1. The inset reveals that the agreement is satisfactory also on a shorter time 
scale. The QUAPI-parameters are K = 4 and At = 0.75 for the symmetric, and At — 0.3 for the asymmetric case, 
respectively. Also for a higher temperature T = 0.2 the agreement is very good, see Fig. ^b.). 

All results exhibit a single exponential decay at long times. This reveals that the bath parameters have been chosen 
such that the dynamics is indeed incoherent and no quantum coherent oscillations can be observed. In absence of 
a static asymmetry e = 0, the (averaged) asymptotic population of the left well is clearly 0.5. This holds for the 
undriven (s = 0) as well as for the driven case (resonant driving s = 1.0,57 = 0.815). However, in presence of a 
bias e = 0.08, the (averaged) asymptotic population of the left well falls below 0.5. The effect of the additional 
time-dependent driving is to increase the asymptotic population on the left, see Figs. ^ and ^ b.). The quantum 
relaxation rate and the asymptotic population of the left well are studied in greater detail in the subsequent Section 
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VI. THE QUANTUM RELAXATION RATE 



The generalized master equation ([Sq ) is an integro-differential equation that governs the decay of the population out 
of one (metastable) welL However, to extract analytically one single rate, which rules the interesting dynamics on the 
largest time scale, requires further approximations. Motivated by the numerical fact that the decay of the population 
is observed to be exponential with one single exponent (see Section VC), we proceed by invoking a Markovian 
approximation for the GME (|5^). This approximation yields a set of coupled ordinary first-order differential equations. 
In absence of external driving, the corresponding coefficients are time- independent. For a driven system, they depend 
on the actual time variable t. When the frequency of the periodic external driving is of the order of the frequency 
associated with the interdoublet energy gap or larger, O > uJo, the averaging of the dynamics over a full driving period 
is appropriate. After averaging, the coefficients of the set of coupled first-order differential equations then assume 
time-independent values. They form the (time-averaged) rate matrix. The smallest real part of the eigenvalues of this 
rate matrix yields the relevant rate, the quantum relaxation rate, which rules the dynamics on the largest time scale. 

It is this novel expression for the quantum relaxation rate which constitutes a second major result of this work. 
On one side, we consider shallow barriers AU > Huiq as well as high barriers AU ^ hojo. In the latter case, the 
condition for the validity of a semiclassical treatment is met. Put differently, since we deal in our approach with 
discrete energy eigenvalues, the semiclassical limit is reached when the number of levels below the barrier becomes 
large. In this case, however, the numerical solution of the GME becomes intractable. On the other hand, we may 
consider temperatures fceT' ~ hcoo, such that the higher lying energy doublets cannot be neglected, as well as lower 
temperatures k^T <^ Hujq. In fact, our analysis contains the spin-boson solution, being the appropriate limit when 
<C hioo and in the absence of strong resonant driving. Finally, we can allow for large driving amplitudes and 
interdoublet resonant driving frequencies. In this latter case, both the restriction to a two-level system as well as an 
equilibrium semiclassical analy sis is prohibited. 

In the following Section VIA, we describe the Markov approximation for the generalized master equation. In Section 
VI B, the quantum relaxation rate is determined as the smallest real part of the eigenvalues of the rate matrix. 



A. Markovian approximation 

The starting point is the generalized master equation (|58|). The inhomogeneities I^(t,to) on the r.h.s. do not 
contribute to the long-time dynamics since they decay exponentially with time on a rather short time scale, see Eq. 
( |63| ) for the inhomogeneity in the case of the double-doublet system determined within lowest order in Aj. Hence, 
this term can be neglected. 

We assume furthermore that the characteristic memory time T,„cm of the kernels of Eq. (|5^) is the smallest time 
scale of the problem (Markovian limit). This means that we can substitute the argument of pp^{t') under the integral 
by the time t, and draw p^vit) in front of the integral. Moreover, the upper limit t of the integral can be replaced by 
oo. We then obtain the Markovian approximated generalized master equation 

M 

P^^^^{t) = E r^,(t)p,,(t) (66) 
i/=i 

with the time-dependent rate coefficients 

/•oo 

r^.(i) = / dTn^,it,t-T). (67) 
Jo 

The explicit time-dependence of the rate coefficients reflects the explicit time-dependent external forcing. In the case 
without external driving, the rate coefficients in Eq. (|67|) become time-independent. 



1. Analytic result for the case without driving 



To obtain specific results, we investigate the lowest order for the kernels Ti^^i/. The time independent rate coefficients 
then read, to lowest second order, 

A 2 /.CO 

r|fj = ^Jdr exp {-{q^ - q^f S{t)] cos [(F, ~ F^)t - [q^ - q,fR{T)\ , 

p^v. (68) 
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The used quantities have been introduced in Eqs. ([ij), (|12|), (Eq), (M) and (B0|). The conservation of probabihty 



requires for the diagonal elements of the second order rate coefficients that v's^J — — J^k^l- ^'f^'^ ■ integral in 

Eq. (^8|) can be solved numerically by standard integration routines | |6^ . However, an analytical soluti on can also 
be derived. In the limit ojr.t ^ oo H, the correlation functions S{t) and R{t) assume the form in Eq. (B3) for the 



(2) 



real part S{t), and in Eq. (B6) for the imaginary part R{t). After some basic algebra, we obtain for the Markovian 
approximated rate coefficients the expression 



p(2) ^ ^A"^ 



hp 



hpujc 

2n 



(69) 



with r(z) being the F-function 



2. High-frequency-driving 

To extract an average long-time relaxation rate in the case with driving, we choose the external driving frequency 
n to be of the order of the intcrdoublet level spacing uJo « wo- This assumption is met, for instance, if one wishes 
to pump population from the lower to the upper doublet by an intcrdoublet resonant field. A time-average of the 
time-dependent rates (Krylov-Bogoliubov-scheme) in Eq. ( |67| ) over the driving period is then appropriate. In general, 
this averaging procedure is reasonable when the driving frequency is much larger than the time scales related to 
tunnehng, i.e., when 3> Af , Af , ... where the Af are the tunnehng splittings of the doublets. 

We insert the explicit shape of the periodic driving sit) — ssin(f7i) in the second order kernels in Eq. (|64|). The 

averaging with respect to the driving frequency reads (yp,vif))a = (f^/27r) dtV The integration over t 

can be performed if one represents the sine- and cosine- function in terms of Bessel functions J„(x) [ p3[ . The only 
remaining non-zero part is the one which contains the zeroth Bessel function Jq(x). The time- independent averaged 
Markovian rate matrix elements to second order emerge as 

rp(2) ^ (rW(t))n 

" io ^^'^^p{^('?a'-9!^)^'5'(t)} Jo f— ((7p-g^)sinf-Tjj 

X cos [(F, - F^,)t - [q^ - q,fR{T)\ , il^v. (70) 

Like in the non-Markovian case, the conservation of probability implies for the diagonal matrix elements the condition 
r^^'^^^ = — X^k/;/ TkI,'^^'. This expression reveals that the infiuence of driving is different for each pair of DVR-states 
since the explicit distance — enters in the argument of the Bessel functions. The averaged rate matrix elements 
cannot be calculated in closed analytical form as in the undriven case; however, they can be obtained numerically by 
standard integration routines . 



B. The quantum relcixation rate 



Since the diagonal elements Pp^(i) obey Eq. (|6^), the long-time dynamics in this regime is ruled by a single 
exponential decay. In the case without driving, the rate matrix F^^ is already time-independent, and cquivalcntly for 
the case of high-frequency driving after the averaging procedure. Both cases reduce to a structure 

M 

PA^MW-Er|.TW(i), (71) 

where the superscript (av) means that the formula holds for the averaged as well as for the time-independent case. 
This set of coupled ordinary first-order differential equations can be decoupled via a diagonalization procedure. If one 
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denotes the elements of the transformation matrix by S^^ and the eigenvalues of the (averaged) rate matrix by A^, 
the diagonalized (averaged) rate matrix reads 

'^)^1K,1^^KIK2^K,1V — ^iJ.^llV ■ (72) 

The general solution of the (averaged) Markov approximated GME is obtained to be 

M 

Pmm(0= V('5"')-e^''^*~*°VKK(io). (73) 

Since r|fj^ is a stochastic matrix, i.e., the diagonal elements of the (averaged) rate matrix are the negative sum of the 
matrix elements of the corresponding columns, one eigenvalue equals zero, i.e., Ai = (conservation of probability). 
Therefore, 

M M 

P,,{t) = + E E S,AS-'Ue''"'''-'"^P..{U) , (74) 



u=2 K = l 



with = J2k=i 'S'/j,i('S'^^)i,kPkk(^o) being the asymptotic population of the DVR-state \qf_i). 

The rate which determines the dynamics on the largest time-scale is the smallest non-zero absolute value of the 
real part of the eigenvalues of the (averaged) rate matrix, i.e.. 



r^'^^) = min{|ReA^|; 1/ = 2, M} . (75) 



It is called the quantum relaxation rate. 



Likewise, the asymptotic population of the left well is readily obtained from Eq. (73). It reads 



L 

i^iSt^E^^M- (76) 

To compare the predictions of the Markovian approximated master equation ( |66| ) with the res ults o f the generalized 
master equation (|5^), we recall the outcomes presented in Figs. ^ and ^ of the previous Section V_C. The Markovian 



results are indicated by the dashed lines. In all investigated parameter combinations, the agreement between the 
generalized master equation, the predictions of the QUAPI algorithm and the Markovian master equation is very 
good apart from minor differences. The rate of the decay is described accurately by r''*^^ as well as the asymptotic 
population of the left well. Also in presence of a time-dependent driving, the averaging yields the correct averaged 
dynamics. This allows for the conclusion that, in the investigated range of parameters, the driven dissipative multi- 
level system is adequately described by the Markovian approximated master equation with the eigenvalues determined 
from second order gNICA. 



VII. RESULTS: QUANTUM RELAXATION RATE AND ASYMPTOTIC POPULATIONS 

In this Section, we present results for the quantum relaxation rate F^**^^ and the asymptotic population P{^^ inside 
the left well for the (driven) double- well potential, Eq. (||). Throughout the following Sections, we choose a set 
of typical dimensionless parameter values. The corresponding dimensionful values follow from the standard scaling 
procedure described in Appendix The barrier height is consistently chosen to be Eb = 1-4. This implies that two 
doublets lie below the energy barrier and the other energy states lie above the barrier, see Fig. ||. Moreover the lower 
tunneling splitting is Af — 3.60 x 10~^, the upper tunneling splitting is Af = 0.121 and the energy gap between the 
two doublets is oJo — 0.815. This choice is mainly motivated by the fact that we explicitly want to investigate the 
intermediate regime between the two-level approximation and the semiclassical regime. Furthermore, such a shallow 
barrier height is convenient for numerical reasons. The splitting of the lowest doublet decreases exponentially with 
increasing barrier height. 

In the following subsections, we investigate on the one hand the double-doublet system, i.e., M = 4, for the barrier 
height Eb = 1-4. We expect that the results are qualitatively similar for larger barrier heights when more than two 
doublets lie below the barrier because the spectrum is then similar to the double-doublet case. On the other hand, we 
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study the question of convergence with increasing number M of energy levels for the case of = 1.4. For Af — > oo, 
the multi-level systems is equivalent to the spatially continuous potential. 

We consider two typical situations: The unperturbed energy spectrum in the symmetric case (e = 0) exhibits avoided 
level crossings, see Fig. |^ a.) below. The second case refers to a tilted potential with e = 0.08 where the energy levels 
are rather strongly separated, see Fig. ^ a.) below. 

The parameters for the time-dependent driving are typically chosen in such a way that the regimes of a weak 
(s=0.1) and a strong (s — 1.0) driving amplitude are covered both. In the first case, the potential stays permanently 
bistable while in the second case the potential assumes intermediate monostable configurations. With E-q = 1.4, the 
critical amplitude where bistability vanishes is Scrit = 0.64. The driving frequency is typically chosen either to be in 
resonance with the interdoublet energy gap, i.e., =TDq = 0.815, or off resonance, i.e., — 0.2. 

The typical choice for the temperature is T = 0.1, being a low to intermediate temperature. We note that the 
semiclassical expression for the cross-over temperature ||^ yields Tco = 0.12, keeping in mind, however, that our choice 
of the barrier height does not obey the semiclassical condition. 

We use an Ohmic spectral density with an exponential cut-off, see Eq. (|^). The damping constant is chosen to be 
7 = 0.1; it represents an intermediate damping strength. We note that the dependence of the results on the damping 
strength and temperature is exponential. The cut-off frequency is always fixed to be Uc = 10.0. 

Finally, we note that in the following symbols such as • and □ are used to label individual plots in the particular 
figures. Their number is not related to the number of calculated data points, the latter being much larger. 

A. Absence of external driving 

We start with the simplest case of the undriven symmetric double-well potential (s = 0, e = 0) and consider the 
dependence of the quantum relaxation rate on the number M of energy eigenstates. M = 4 denotes the double-doublet 
system. Fig. ^ shows the result for different damping constants 7. A convergence of the rate F can be observed for an 
intermediate damping strength 7 = 0.1 (*). We recall that the lowest tunneling splitting is two orders of magnitude 
smaller than 7 and that the upper tunneling splitting is of the same order of magnitude as 7. For a larger damping 
constant 7 = 0.5 (□), convergence is also obtained. However, for the case of very strong damping 7 = 1.0 (A), the 
result for M = 8 does not coincide with that for M — Q. This fact is due to the following feature: As it can be seen 
from Eq. (^) for the second order kernels of the generalized master equation, the damping constant 7, which enters 
via S{t) and R{t), is multiplied by {q^ — q^Y being the square of the tunneling distance between the two involved 
DVR-states. Upon increasing the number M of energy levels, the DVR-eigenvalues lie more dense in position space. 
Hence, some distances become small and the multilevel system effectively flows to weak damping. The small effective 
damping is no longer sufficient to suppress long intervals in the off-diagonal states. Thus, the gNICA in second order 
is no longer applicable and contributions of higher orders of the integral kernels in the GME have to be taken into 
account. A more detailed discussion of the effect of the flow to weak damping is postponed to the Appendix ^ 

The question of convergence of the quantum relaxation rate with increasing M in presence of time-dependent driving 
is investigated in the sections below. 

B. The influence of external (time-dependent) driving forces 

In this Subsection we investigate the role of external driving. This external perturbation can be either a static 
potential asymmetry (bias) or a time-dependent periodic driving (ac-driving), or simultaneously both parts are 
present. 

1. Dependence on a static bias, no ac-driving 

Adding a static asymmetry renders one (in our case the left) of the two formerly stable potential minima a metastable 
minimum. The consequences for the spectrum of the bare system are that avoided level-crossings occur for particular 
values of the asymmetry, see Fig. ^ a.). At such avoided level-crossings tunneling is usually enhanced. This effect is 
known as resonant tunneling. This situation, however, is modified in the presence of a moderate to strong damping. 
The case of a strong system-bath coupling is considered in Fig. ^ b.) where the relaxation rate shows peaks at 
particular values of the static bias. Their position strongly depends on temperature. At low temperatures T = 0.05 
to r = 0.15 (full lines), we observe three relative maxima at e = and around e = 0.12 and e = 0.25. First, we 
emphasize that the quantum relaxation rate initially decreases when the bias is increased from zero, i.e., when the 
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effective barrier height is decreased. This feature is a typical quantum mechanical footprint. In a classical system, the 
relaxation rate grows when the barrier is lowered Second, we note that the two peaks at nonzero bias values are 
shifted to smaller values of the asymmetry compared to those bias strengths where the avoided level crossings occur, 
see Fig. ^a.). This indicates that we are no longer in a weak-coupling regime but encounter already strong incoherent 
tunneling. 

Increasing the temperature results in a decrease of the amplitude of the peak at e = 0.25. This indicates the 
enhanced destruction of the resonant tunneling phenomenon. Moreover, the peaks broaden. At an intermediate 
temperature T = 0.3 (dashed line) one characteristic peak occurs at e = 0.2. Its height is smaller than in the low 
temperature cases which indicates that tunneling is reduced compared to the low temperature case. This is mainly 
due to the enhanced environmental level broadening. However, we still observe a clear decrease of F when the bias is 
increased from zero onwards; therefore we conclude that quantum tunneling still occurs. This, however, is no longer 
observable for the high temperature case T = 0.5 (dotted line). A relaxation rate which grows with increasing bias is 
a signature of classical behavior. 

The question of the convergence of the rate with an increasing number M of energy states is addressed with Fig. 
0. We show four different cases for M = 6 (•) and M = 8 (□) and T = 0.1 (fuU line) and T = 0.15 (dashed 
line), respectively. The low temperature case T = 0.1 shows a clear convergence when M is increased from M — 6 
to Af = 8 for small asymmetries up to e = 0.15. For larger asymmetries, the two results, however, do not agree. 
This behavior can be resolved as follows: For large asymmetries, the left well is strongly lifted above the right well. 
Moreover, the position eigenvalues on the left side move towards each other and are densely located. This means that 
the tunneling distance of the corresponding transition becomes smaller which in turn reduces the effective damping. 
Then a flow to weak damping occurs and the second order gNICA breaks down. The same explanation holds for the 
larger temperature T = 0.15 where the results for M = 6 and M = 8 show qualitatively a similar behavior up to 
e = 0.15. 



We investigate in the following the asymptotic population P{^^ of the left well determined in Eq. (7^). For the 
symmetric case, it assumes the value 1/2. In presence of a positive asymmetry e > 0, P{^^ is smaller than 1/2 
since the left well is energetically higher. Fig. | shows P{^^ as a function of e for two different temperatures (solid 
line) for M — 8. The damping constant is chosen to be 7 = 0.1. For comparison we additionally show the asymp- 
totic population obtained from a Boltzmann equilibrium distribution for the same parameters (dashed line), i.e., 
p(oo) = exp(— Ho/fcer). First, we note that the asymptotic population decreases exponentially with increasing bias. 
Second, we emphasize that the often made assumption of a Boltzmann equilibrium distribution is only valid for an 
infinitesimally small coupling of the system to the environment . The depicted results for M = 8 have already 
been converged and are not distinguishable on our scale from the Af = 6 case. This is indicated in the inset of Fig. 
^ for a fixed asymmetry e = 0.08 for two different temperatures. Convergence is also found for the entire considered 
parameter range of e (not shown). The full line again shows the result obtained from gNICA in second order while 
the dashed line marks the results from a Boltzmann equilibrium distribution. 



2. Dependence on the static bias in presence of external ac-driving 

The influence of a time-dependent periodic driving on the quantum relaxation rate is elucidated with Figs. ^ - 
For the case of off-resonant driving. Fig. |^ exhibits a non-monotonic dependence of the averaged relaxation rate 

F'^'^ on the bias. For increasing Af the results approach each other. However, a complete convergence as observed in 

the undriven case (Fig. 0) is not obtained. 

Tuning the driving frequency Vl into resonance, the results in Fig. |l^ show for a fixed driving strength a characteristic 

peak, being almost independent of the temperature T . The position of the peak is sensitive to the driving strength. 

This indicates that the population of the upper doublet is mainly the result of driving and not due to thermal 

population. 

Furthermore, we draw the reader's attention to the strongly (s = 1.0) driven symmetric case e = 0. The low 
temperature relaxation rate for T = 0.1 is larger than for the two other cases with higher temperatures. This is 
opposite to the situation without driving, see Fig. |^ b.) where F is smaller for T = 0.1 compared to T = 0.5. This is a 
typical footprint of a quantum effect: The resonant driving (f2 — 0.815) in the symmetric potential transfers population 
to the upper doublet where tunneling is enhanced because the tunneling splitting is large and the temperature is not 
high enough to destroy coherence completely. 

The problem of convergence of the results with increasing Af is addressed in Fig. |ll| for the resonantly driven case. 
Shown are four different combinations for Af = 6 (•) and Af = 8 (□) for two temperatures T = 0.1 (full line) and 
T — 0.15 (dashed line). We ascertain that no convergence is obtained upon increasing Af . The difference between the 
results for Af = 6 and Af = 8 for this large driving frequency is larger than in the case of off-resonant driving, see Fig. 
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g. Obviously, more than only a few energy eigenstates are necessary to describe the resonantly driven double-well 
potential accurately. This is not astonishing, however, because the driving frequency is in resonance = 0.815) and 
many higher energy levels are excited. Moreover, we stress that the results vary within the same order of magnitude. 

The asymptotic population Pj^^ of the left well is shown in Fig. ^ as a function of the static bias e for three 
different temperatures. The driving frequency is in resonance with the interdoublet energy gap. Compared to the 
strict monotonic behavior of the undriven case, see Fig. H, the results with strong driving (s — 1.0) show a non- 
monotonic dependence on the static bias with a distinct maximum at e = 0.15 being nearly independent of T. This 
maximum has a value of P{^^ ~ 0.5 indicating an equal population of both the metastable as well as the stable well. 
The position of this local maximum is invariant under the choice of the driving strength, as indicated by the weak 
driving result for (s = 0.1). We note that a net population inversion can be achieved for the parameter combinations 
M — 4:,T = 0.1, s — 1.0, £ — 0.14. It is interesting to see that convergence of Pj^^ with increasing M is obtained for 
asymmetries e > 0.2, see inset of Fig. For smaller values of the bias the qualitative behavior in the case of M = 6 
is similar to the case of M = 8. However, the local maximum around e = 0.15 in the double-doublet case M = 4 
vanishes upon increasing M. 



3. Dependence on the driving strength 

Fig. ^ depicts the averaged rate as a function of the amplitude s of the ac-driving field with resonant driving 
frequency fl —ZJo ^ 0.815 for the symmetric double-doublet system (e = 0, AI — 4). Shown are the results for three 
different temperatures. The asterisks * mark the results of an exponential fit to QUAPI results (not shown) and 
confirm the validity of our new analytical approach. 

The averaged rate for the case of a high temperature T = 0.5 is reduced compared to the undriven situation where 
pav^^ r) has a maximum. Upon decreasing the temperature, the relative maximum at s « 0.9 grows out to a glo bal 



maximum for T = 0.1. This resonance is useful for practical applications ("Hydrogen subway", see Section I A) if 
one desires to accelerate the transfer of population from the left to the right well. So not only a resonant driving 
frequency O = ZUg but also a suitably chosen driving strength is important to maximize the transfer. The behavior of 
F'^'^ vs. the driving amplitude is shown in Fig. |lj for an increasing number M of energy states (note the logarithmic 
scale). For small driving strengths (up to s « 0.2), the result for M — 10 does not significantly differ from the case 
for M = 8, indicating numerical convergence. For intermediate to strong driving, the differences increase. However, 
we stress that the results for M — 8 and M = 10 remain within the same order of magnitude. 



4- Dependence on the driving frequency 

The dependence of the averaged relaxation rate on the driving frequency is shown in Fig. |l^ for the symmetric and 
the asymmetric double-doublet system M — A. The results can be viewed as a scan of the spectrum of the driven 
dissipative double-doublet system. At some values of fl the transition from the left to the right well is enhanced. 
The intermediate damping constant 7 = 0.1 leads to a considerable broadening of the energy levels involved in the 
transitions, as can be deduced from the rather broad resonance lines. The symmetric (asymmetric) case reveals a 
distinct peak ai fl = 0.4 {il = 0.5) together with sidebands at the corresponding fractions of fl. The behavior of F'*^ 
for an increasing number M of energy levels is depicted in Fig. |l^ a.) for the symmetric case and in Fig. |l^ b.) for 
the asymmetric case. It can be seen that the additional energy levels yield additional resonance lines. However, if one 
chooses the driving frequency sufficiently far from any resonance line, convergence can be achieved. 

The asymptotic population Pj^^ of the left well as a function of the driving frequency H. is depicted in Fig. ^ for 
an increasing number M of states. Clearly, the results do in general not converge with growing M . This confirms 
our result that for an accurate description of a strongly driven quantum system more than only a few basis states are 
required. 



C. Dependence on the bath parameters 



1. Influence of temperature 

The dependence of the quantum relaxation rate on temperature is depicted in Fig. ^ for the case of the double- 
doublet system M = 4. Shown arc four cases without ac-driving (s = 0), with resonant ac-driving (s = 0.1, il — 0.815), 
without bias (e — 0) and with bias (e = 0.08). We first concentrate on the undriven case s = (full line). Interestingly 
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enough, we find in the low temperature regime (see inset of Fig. ^8|) that the rate first decreases when the temperature 
is increased in presence of a static bias. This is a characteristic quantum feature: In contrast to a classical behavior, 
the quantum relaxation rate first decreases with increasing temperature due to an enhancement of decoherence. Then, 
however, the rate starts to increase again as soon as the higher doublet becomes thermally populated. This typical 
behavior has also been observed experimentally in the context of tunneling of impurities in solids |3^-^ (see also 
Section lA and especially Ref. |39| ) . For the intermediate temperature regime the comparison between the symmetric 
and the asymmetric case reveals another interesting characteristics: One could argue that the almost linear increase of 
the rate with temperature reveals a classical Arrhenius behavior. This, however, is not the case. We are still in a deep 
quantum regime, since F in the asymmetric case is smaller than in the symmetric case! This is again a clear sign of 
quantum mechanics since the symmetric potential with e = corresponds to a resonant tunneling situation. A finite 
bias of s = 0.08 implies a non-resonant tunneling situation. There, the transfer via resonant tunneling is suppressed 
and the rate becomes smaller. In the regime of intermediate to high temperatures, this behavior is inverted, i.e., the 
rate for the asymmetric case with reduced barrier height is now larger as compared to the symmetric case. 

In presence of a weak resonant ac-driving s = Q.l,fl = 0.815, the quantum relaxation rate F^^ for the symmetric 
case £ = at low temperature is larger than the corresponding undriven relaxation rate. Coherent excitations to the 
upper doublet where tunneling is enhanced dominate at low temperature. The presence of an additional static bias 
e = 0.08 renders the averaged relaxation rate F^'^ almost independent of temperature. 

The question of convergence of the rate with increasing the number M of energy eigenstates is addressed in Fig. 
[l9| a.) for the undriven case s = and in Fig. |l^ b.) for the off-resonantly driven case s = 0.1, SI = 0.2. In absence of 
ac-driving, a satisfactory convergence between the case M = 6 and M = 8 is achieved in the low temperature regime 
for both the symmetric and the asymmetric potential. Clearly for higher temperatures the agreement is worse because 
now the higher lying energy states are not negligible. Also in presence of an off-resonant ac-driving, convergence is 
achieved at low temperatures. 

The asymptotic population P{^^ of the left well as a function of the temperature T is shown in Fig. ^ for an 
increasing number M of states. The results for the undriven case s = reveal a satisfactory convergence over the 
entire temperature regime. For comparison, we also depict the asymptotic popul ation ste mming from the assumption 



of a Boltzmann equilibrium distribution. A similar argumentation as in Section VII B 1 (see Fig. |8|) holds to explain 
the disagreement. 



2. Influence of damping 

Fig. Ell depicts the (averaged) rate vs. 7 for an increasing number AI of states. We find for the undriven case s — 
in Fig.^l] a.) convergence for 7 > 0.08 for both the symmetric and the asymmetric potential. However, for smaller 7 
notable differences occur which indicate that the gNICA to second order is not reliable because the effective damping 
becomes too weak (effect of flow to weak damping, see discussion in Appendix ^). The situation is similar in presence 



of an off-resonant ac-driving s = 0.1, fl = 0.2, see Fig. 21 b.). The results for the case of a resonant driving Q — uja 
are qualitatively similar (not shown). 



VIII. CONCLUSIONS AND OUTLOOK 



A novel scheme to investigate analytically as well as numerically tunneling and vibrational relaxation in a strongly 
driven bistable potential was presented. A necessary first step in our approach is the reduction of the system dynamics 
to the Hilbert space spanned by the M lowest energy eigenstates of the static bistable potential. Because the coupling 
to the heat bath is bilinear in the system and bath coordinates, the convenient basis to perform calculations consists of 
the eigenbasis of the position operator, i.e., the so-termed discrete variable representation (DVR). It is this DVR basis 
which permitted us to derive a set of non-Markovian generalized master equations (GME) for the diagonal elements 
of the reduced density matrix. In the studied regime of temperature and damping, the Markovian approximation to 
the GME yields novel analytical results. They agree well, both with those of the full GME and of precise ab-initio 
numerical path integral calculations. In turn, the quantum relaxation rate could be extracted from the Markovian 
rate matrix. The dependence of the quantum relaxation on the five most relevant model parameters, namely bias 
strength e, driving strength s, driving frequency 17, temperature T and damping 7, was outlined in detail. 

We have identified several quantum mechanical footprints in this strongly damped system. The four most pro- 
nounced quantum features are: 

(i) In absence of ac-driving we find resonant incoherent tunneling. This is demonstrated by striking resonances in 
the relaxation rate at distinct values of the dc-bias. 
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(ii) We observe a decrease of the relaxation rate as the effective barrier is lowered by a static bias. This finding is due 
to a reduction of tunneling since the energy gaps forming the tunneling doublets increase with increasing asymmetry. 
In contrast, the relaxation rate in a classical system always increases for a reduced barrier. 

(iii) The ac-driven quantum system shows distinct resonances of the relaxation rate for particular values of the 
driving amplitude. Especially at low temperature, the relaxation rate is enhanced by driving as compared to the 
undriven case. 

(iv) A non-monotonic dependence of the relaxation rate on temperature is observed. Increasing the temperature 
in a classical system always increases the rate. However, in a quantum system, a higher temperature induces a larger 
population of the energetically higher lying doublets, where tunneling is favored. Increasing temperature further 
renders the quantum system more incoherent and the relaxation via tunneling is again hampered. The rate therefore 
decreases before it grows again due to thermal hopping. 

Our analysis furthermore permits to determine the asymptotic population of the left metastable well. We have 
shown explicitly that a Boltzmann equilibrium distribution (both in absence and in presence of an ac-field) is not 
attained for the chosen set of parameters. 

The GME in Eq. (^sj) and its time-inhomogeneous Markov approximation in Eq. (^6|) treat the external driving 
field exactly, and are reliable for moderate temperatures and/or moderate damping strengths. Indeed, the equations 
become exact for Ohmic dissipation at high temperatures. Thus, our approach complements a Redfield-type analysis 
being appropriate for weak damping. In contrast to semiclassical calculations we can consider shallow potential 
barriers substaining only a few doublets below the barrier. Moreover, in contrast to semiclassical imaginary-time rate 
calculations, we are not limited by the requirement of thermal equilibrium at adiabatically varying external fields. 

A major restriction of the presented method is that the generalized non-interacting cluster approximation (gNICA) 
that we used to obtain the GME turns out to be useful in praxi for numerical purposes only when (i) the number M 
of levels remains moderately small (M < 10), and when (ii) the truncation of the gNICA kernels in Eq. ( ^9| ) to lowest 
order in A^^ is appropriate, see Eq. (^^. Clearly, the number M of levels involved in the dynamics increases with 
increasing temperature and/or for large driving strengths s, and/or for resonant driving frequencies ft. 

Taking into account a larger number Af of basis states implies that the position eigenvalues lie more dense in position 
space (in the limit of infinitely many energy eigenstates the distance between neighboring DVR-points is infinitesimally 
small). However, we have observed in the preceding Sections that the square of the distance between two DVR-points 
enters as a prefactor for the twice integrated bath autocorrelation function S{t) + iR{t) in the second order kernels 
in Eq. (64) of the generalized master equation. This implies that upon increasing M the effective damping of each 
transition is reduced and the multi- level system effectively flows to weak damping. For small effective damping the 
noise action does no longer suppress long intervals in the off-diagonal states (clusters) and higher than second order 
transitions start to contribute. 

To deal with this effective weak coupling situation (which occurs for large M even when the global coupling constant 
7 = r]/Ai is not small), a procedure similar to the one used by Zwerger |^ and Grabert, Ingold and Paul to 
investigate transport in Josephson junctions can be used, see Appendix ^. 

Due to its very general nature, the newly developed analytical technique contains a large potential for applications 
to speci fic experimental situations. Several possible applications for experimental systems have been discussed in 



Section [A 



A prominent question for future work refers to the behavior of the crossover to the classical regime. Our results 
should merge into those of the quantum Kramers rate ||] for semiclassical barrier heights, i.e., AU/{huJo) = Eb ^ 1, 
in the case without ac-driving. Once this regime is explored, the method can be generalized to time-dependent 
semiclassical quantum systems. However, one has to be aware that for the driven system a chaotic dynamics generally 
occur |6^. We expect that the dynamics in the semiclassical regime involves an increasing number M of DVR states. 
Then the effect of the flow to weak damping occurs. The analysis in Appendix ^ represents the starting point for 
future investigations towards this challenge. 
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APPENDIX A: SCALING TO DIMENSIONLESS QUANTITIES 



For the specific calculations, we introduce in this Appendix dimensionless quantities. They are obtained by scaling 
the Hamiltonian in Eqs. (]l|) and and the environmental parameters specified in the Eqs. (|]) and (^). The relations 
read 



i = iuot, q{i) = ^Mujn/hq{t= Eb = AU/hwn., 



We omit all the tildes for the sake of better readability. 
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APPENDIX B: THE BATH CORRELATION FUNCTION 

In this Appendix we give the explicit expressions for the twice integrated bath correlation function Q{t) = S{t) + 
iR{t) of Eq. ( |Tl| ) with an Ohmic spectral density of Eq. (||), see also The real part S(t) can be evaluated analytically 
by solving the integral in terms of the ■^-function. One arrives after some algebra at the exact expression 



S{t) 



In 
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For finite temperatures k^T <^ hwc {scaling limit), we obtain with r(l + i?;)r(l — iy) — Try/ sinh(7rj/) the expression 
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This illustrates that the correlations between the paths are damped out exponentially at long times for a low tem- 
perature Ohmic bath. The temperature-independent imaginary part R{t) in Eq. ( [Tl| ) can be determined exactly. We 
obtain after the integration over lu 



R{t) — — arctan(a;ci) ■ 

TT 

For long times uict — > c», the arctan-function approaches the Heaviside function, i.e., 

R^t) "^^^^ |e(0 , 

so that the imaginary part R{t) becomes a constant function for all times t. 
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APPENDIX C: NUMERICAL ITERATION SCHEME FOR SOLVING THE GENERALIZED MASTER 

EQUATION 



The generalized master equation ( |58| ) is a set of M coupled integro- differential equations with inhomogeneities . For 
the one-dimensional case M = 1, the GME would be of Volterra type. For this case, several standard techniques for 
a numerical treatment are known |p3,c2| and common numerical libraries [p2| supply codes for this case. The general 
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M-dimensional case is far from being standard and we are not aware of available algorithms in the literature. The 
non-locality in time detains us from diagonalizing the kernel rate matrix and thereby decoupling the equations. We 
develop in this Appendix a rather simple numerical algorithm for solving the general set of M coupled inhomogeneous 
integro-differential equations. 

We start by formally integrating the GME (|5^) once and choose the integration constants such that the initial 
conditions are fulfilled. After interchanging the order of integration, we obtain 

M pt nt 

Pf.f.{t)^Yl / dt' JC^,{tX)Puu{t')+ dt'I^{t',to) + p^^{to), (CI) 
where we have defined the integrated kernels 

t 



jc^,[tX)= / dt"n^,it",t'). (C2) 



M pt+At pt+At 

,(t + Ai) = E / dt'ICf,,{t + At,t')p,,{t')+ dt' I^{t',tn)+Pf,f,{to) 



In the next step, we iterate the integrated GME from time t to time t + At, and split the integrals to obtain 

M pt+At pt+At 

Pf^f^C ' ' " 

Y dt' dt"nf,,{t",t')p,,{t')+ / dt' If,{t',to)+p^,^{to) 

y=X J to Jt' J to 

M pt i-t+At 

dt"H^,{t",t')p,,{t') 



i-t i-t 

u=l "'to 



/t+At M pt+At 

dt' I,,{t',to) + Y^ dt' JC^,{t + At,t')p,,{t') 



v=l "^t 

M pt nt + At 



= pmm(^) + E / dt'p,,{t') dt"n^,{t",t') 

y=l Jt 

pt+At M pt+At 

+ dt' I,,{t',to) + Y, dt' lC^,{t + At,t')p,,{t'). (C3) 

So far, every manipulation was an exact transformation. To proceed, we have to invoke an approximation for the last 
term in Eq. (|C^), i.e., that one involving 1C^y{t + At, t'). It is this term only which requires the knowledge of PuA^') 
in the time interval \t,t + At]. All the other terms only need Pvv{t') up to time t which are known. For that reason, 
we use in a third step the simplest approximation rule for the integral, i.e., the Simpson trapezoid rule, to obtain 

M pt+At 



V / dt'ICf,^it + At,t')p^^{t')!i^ 
E ^{^/^-(* + At, t) p^^t) + JC^,{t + At,t + At) p,,{t + At)} . (C4) 



The corrections are of the order of At"^. With Eq. (C2), it follows that /C^^(t + At, t + At) — and we arrive at the 
final iteration scheme 

At, 



(t + At) = E P'^'^ { -^M^ + ^^M^ (t + At, t) } 

M pt+At pt pt+At 

+ E / dt" dt'n^,{t",t')p,,{t')+ dt' i^{t',to).. 



(C5) 



where the is the Kronecker symbol. We note that this iterative procedure requires the knowledge of Pwif!) in 
the time interval t' G [to, t] when propagating from time t to time t + At. Furthermore, we remark that this iterative 
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algorithm is not restricted to the lowest order for the kernels 'Hfj,u{t,t') and the inhomogeneities /^(t,to)- Finally, we 
observe that the integrations from t to t + At for each step can be perform ed n umerically to a very high precision 
such that the only relevant numerical error arises from the splitting in Eq. (C4). Practical calculations reveal that 
the time step At has to be chosen rather small since the problem is similar to a stiff differential equation. In praxi, 
this means a value of the order of At — 10^^ or smaller. This rather small value for At restricts the applicability 
of this very simple and straightforward iteration scheme to problems where the decay is not too slow. More refined 
numerical algorithm are imaginable which could circumvent this problem. 

Because we are interested in the long-time dynamics, iterations up to times t = 5000 can be necessary. Since the 
kernel matrix elements contain exponentials with asymptotically linearly growing exponents (see Eqs. (^l]) and (B4) 
in Appendix^), the memory ranging from time t backwards to time to can be cut-off after some fixed time span Tr- 
Then, the memory is only relevant over the time interval t — and all exponentially small contributions from the 
time interval [to,t — Tr] can be neglected. This accelerates the iteration considerably and avoids too large arrays for 
the storage of the intermediate values of Puiy{t'). However, lowering the temperature demands an increasing memory 
range r^. 

Once the diagonal elements Pfifj,{t) arc known, the population of the left well Picft(i) can be evaluated according to 
Eq. (H). 



APPENDIX D: EXAMPLE: A SINGLE PATH SUBJECT TO DISSIPATION 

The purpose of this Appendix is to illustrate the general derivation of the cluster function, Eq. (^), within the 
gNICA scheme introduced in Section pA^ . This approximation scheme is the basis for the derivation of the generalized 
master equation presented in Section^ For simplicity, we pick one single path subject to dissipation and starting in 
the diagonal state (/xq, po) = (qk, qj. = qk) and ending after N = 5 jumps in the diagonal state (/X5, /15) = (q„, = 
The full path is illustrated in Fig. p2| a.) in a time-resolved picture and in Fig. b.) as a path jumping between the 
states of the reduced density matrix in the (g, q')-plane. The path is characterized by the sequence of index pairs 

(MO,/^o) {Pl,l^l) (A*2,M2) {^,'^3) (M4,^'4) (^^5,^5) /p-^N 

= iik.q'i) {quqi) iqui'm) iihin) {qn,qn) ■ 

It contains three sojourns (from to to ti, from ^2 to ^3 and from ts to t), cf. Fig. a.). Moreover, the path contains 
two clusters (from ti to ^2 and from ^3 to ^5). The details of the visited states are illustrated in Fig. ^ b.): The 
diagonal states are marked by diamonds O and the visited off-diagonal states by filled circles •. For our purpose here, 
it is only important to distinguish between diagonal and off-diagonal states. The specific indices of the states are 
irrelevant. 

We evaluate now the contribution X{t) of this specific path to the full path sum in Eq. (|3^). The product of 
the factors Aj in Eq. ( |3^ ) yields a proportionality factor which we omit for the moment for simplicity. For further 
convenience, we consider the time-mdependent system, implying that the diagonal elements E^^., see Eq. (|30|), in the 
Hamiltonian Hg^^ are constant in time. The generalization to time-dependent systems is discussed in Section VA. 
Therefore, we introduce the short-hand notation in Eq. ( |33| ) according to 

AEj = Efj,^ - Ey^ , 

Ti,{ti - tj) = c^p{CiS{ti - tj)Cj + i^iRiti - t,)xj} . (D2) 



The contribution of the specific path to the path sum follows from Eq. (33) as 

rt pta Hi nts nt2 

I{t)= / dts / dti / dts / dt2 / dti eyiY>{i[AEi{t2 - ti) + AE:i{ti ~ t:i) 

J to J to J to J to J to 

+AE4t5 - t4)]}^i,o(ii - io)^2,o(i2 - ^0)^24(^2 - ^1)^3,0(^3 - io) 

X-^3,l(^3 — tl)^3,2{t3 — t2)^4,o{t4 ~ to)-^44(i4 — ^l)-^4,2(i4 — ^2) 

X-^4,3(i4 — i3)-^5,o(i5 " io)-^5,l(*5 " ^l)-^5,2(i5 " i2)-^5,3(i5 " t3)^5,4{t5 — ^4) • 

(D3) 

Eq. ( |D3| ) is still exact. We apply then the generalized non-interacting cluster approximation gNICA in different 
steps. 
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(i) Neglect of the intercluster correlations. 
Let us consider the product 



Vi = T^sih - - t2)Tis{U ~ ti)Tia{U - t2)^5a(i5 - ^1)^5,2(^5 - ^2) • 



(D4) 



If we assume a very large cut-off frequency, i.e., Wc ^ 00, we can approximate the real part S{ti — tj) of the bath 
correlation function by its linearized form, Eq. (B4), and the imaginary part R{ti — tj) by the constant value 77/2, see 
Eq. (B6). Then, Vi can easily be brought into the form 



Vi » exp <; [^3 + 6 + 6] 
V 



hi3 



(i3-ti)a + 1^(^3-^2)6 + -In 



hp 

[6+6; 



27r 



(6+6) 



1. 



(D5) 



In the last step, we have used the neutrahty of each cluster, i.e., 6+6=0 and 6+6+6 = 0. The overall result 
is that we can set in the influence functional the product with all the couplings between different clusters equal to 1, 
i.e., we disregard intercluster correlations. 

(ii) Neglect of cluster-sojourn correlations. 
Next, we consider the product 



^2 = •?^3,o(^3 — io)-^4,o(^4 — io)-^5,o(^5 " ^o) — 

exp{^35(t3 - io)6 + US{U - io)6 + 65(^5 - ^0)6 
+i[i?,R{h - to)xo + 6^(i4 - io)xo + 6i?(i5 - io)xo]} 



(D6) 



describing the interaction of the clusters with the initial state characterized by (xo;6)- Since we start in a diagonal 
state, it follows that 6 = 0- Moreover, we apply the same argumentation like in (i) for the imaginary part R{t) and 
obtain 



7^2 « exp {z [^3 + 6 + 6] fxo} = l. 
The corresponding argumentation holds for the third product J-2fi{t2 — io)-?^i,o(ii 



(D7) 



to). 



Step (i) and (ii) contain all the correlations disregarded within gNICA. In contrast, we entirely keep the m<racluster 
interaction J^^^iit^ — ^4) as well as the interactions of the particular clusters with the directly preceding sojourn, i.e., 
•^4,3(^4 — ^3) and ^^^^{1^ — t'i) and J-2,i(t2 — ti)- After reordering the integrals we obtain 



to 



xT^Ah - ^3)^^5,3(^5 - ^3)^^5,4(^5 - u) 



dh 



to 



dti exp [iAEi{t2 - h)] ^^2,1(^2 - ^i) 



(D8) 



to 



This expression can be treated more conveniently after a Laplace transformation to X(A) = Ct{T{t)} = 
dtexp{—Xt)l{t). Using the property Ct{J^^dt^f{t^ — to)} = j£j{/(i)} the integration over t^ yields the 
factor 1/A. The remaining function f{t) to be Laplace transformed can be recast into the convolutive form 
f{t) = Jl dt3gii~t3)h{t3) with g{i-t3) = dt4eMi[^E3{ti^t3) + AE4{i-ti)]}TiAti-t3)T5Ai-t3):F5A{i-ti)} 
and with h{t3) = j^'^ dt2 f^^ dti exp[iAEi{t2 — ti)]J^2,i{t2 ^ ti)- By application of the convolution theorem 
Cj:{J^^ dt^git — t3)h{t3)} — Ci{g{i)} Ci{h{t)} and by performing the integration over t2^ we obtain the product 
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xi^-j^ t^^i exp[iA£;i(t-ti)]^2,i(t-ii)| ■ (D9) 

The Laplace transforms are the contributions of the full intracluster interactions of the two clusters. Performing the 
Laplace transforms and transforming to the time differences Tj = ij+i — tj, we finally arrive at the expression 

JgNiCA(^) = / drj dT3 exp{-A(T3 + t^)} exp{i(Ai?3r3 + AE^n)} 
^ Jo Jo 

X•?^4,3(•r3) •?^5,3(t3 + T4) T^^iir^) 

1 1 

x-y c^Tlexp{-ATl}exp{^A^;lTl}.F2,l(Tl)-, (DIO) 

where the last factor 1/A appears after the integration over the first sojourn ti — <o- Eq. ( pDlOj ) is an example of a 
contribution of one specific path to the total path sum in Eq. (^). 

APPENDIX E: HARMONIC WELL APPROXIMATION 

In this Appendix we give explicit results for an approximation for the energy eigenstates in the wells of the double- 
well potential (||) without external forces, i.e., £ = and s = 0. The scaling o f th is Hamiltonian is performed 
according to the standard procedure described in the previous Appendix ^ see eq. ( |Al[ ) . In the literature , often 
the assumption is made that the energy eigenvalues and the localized states in the two wells can be approximated 
by those of a harmonic oscillator potential whose minimum coincides with the single well minimum. The localized 
states of the double-well potential are linear combinations of the sy mmetr ic and the antisymmetric energy eigenstate 
corresponding to one doublet. They have been introduced in Section [II D| , Eq. for the case of the double-doublet 



system and a generalization to more than two doublets is straightforward. 

The eigenenergies and the energy eigenstates of a spatially shifted harmonic potential are given in dimensionless 
units by 



and 



£n^n + l/2, n = 0,l,..., (El) 



Mq) = {n\q) = (2" exp \ -^{q - qof } Hn{q - go), n = 0, 1, . . . , (E2) 



where go = ±v8^b is the position of the minima of the double- well potential with barrier height Eb and Hn{q) are 
the Hermite polynomials. 

Fig. p3| shows the results of this approximation for two cases of a barrier height E^ in the deep quantum regime, i.e., 
i?B = 1-4 (Fig. a.) with two doublets below the barrier and E-q = 2.5 (Fig. |23|b.) with three doublets below the 
barrier. For comparison the exact (numerically obtained) wave functions are depicted. In both cases, the agreement 
is not convincing, as expected. Increasing the barrier height improves the agreement for the lower lying states. The 
states lying closed to the barrier top show a noticeable disagreement. Especially the part of the wave function which 
penetrates the barrier and which is in turn associated with tunneling is underestimated. 

We note that for graphical reasons, the harmonic states are positioned at the exact eigenenergies of the double-well 
potential. However, in the approximation the harmonic eigenenergies are also shifted compared to the exact one (see 
scaling on the ordinate). To be specific, the interdoublet energy gap in the case of E-q = 1.4 is cJq — 0.815 in contrast 
to ZUq = cjQ = 1 for the corresponding harmonic potential. For the case of E^ = 2.5 the lower interdoublet splitting 
is cJo,i = 0.892 and the upper interdoublet splitting is IUo,2 = 0.805. These values have also to be compared with 
Uq = t^o = 1 for the corresponding harmonic potential. 
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The deviations up to 20% are certainly not astonishing since the condition for the apphcabihty of the harmonic well 
approximation is a rather high barrier, i.e., Eb 3> 1. This is also the relevant condition where semiclassical methods 
are applicable to calculate quantum relaxation rates |^ . They are rather simple compared to existing approaches Q . 
For barrier heights of the order of the interdoublet spacing, this approximative treatment of the eigenenergies and the 
wave functions is not applicable. 



APPENDIX F: FLOW TO WEAK DAMPING 



In this Appendix we describe the main steps of a proposal of how to deal with the effect of the flow to weak 
damping. A deta iled investigation is still in progress. This part is to be viewed as a starting point for future work 
(see Section VIIl) and should only point out that this behavior can also be treated within the formalism of real-time 
path integrals. 

To deal with this effective weak-coupling situation, the twic e int egrated bath autocorrelation function S{t) + iR{t) 
in the asymptotic limit of the scaling limit, i.e., Eqs. (B4) and (B6) of Appendix may be used. The deviations from 
the exact form of the autocorrelation function are small for high temperatures and/or weak damping. We use these 
approximative expressions in the kernels in Eq. ( p9|) for the generalized master equation and obtain for the discrete 
influence phase (see also Eq. (pF 



N l-l N l-l 

1=1 j=0 



1=1 3=0 

N l-l 

— ^EEe'(^'-^.)e.-Jln 

1=1 3=0 



ll^f ^^^^ 
TT V 27r 



■ V 



N l-l 



1=1 j=0 



IX j 



N l-l 

EE^'^. 

1=1 3=0 



(Fl) 



We define 



P3 



1=1 



(F2) 



and note that each cluster has a cumulative path weight of zero, see Eq. 
elementary rearrangement of the double sums in order to obtain with tj — tj 



|), i.e., pn = 0. This allows for an 
tj^i the expressions 



N l-l N 

EE^'(^'-^^)^^^E^^^^?' 

1=1 j=0 3 = 1 

N l-l N-1 

EE^'^j = " E ^iPi- 

1=1 j=0 3 = 1 



N l-l 

EE?'^. 

1 = 1 j=0 



1 ^ 



(F3) 



Inserting these equations into the kernels, Eq. (|59|), it follows that 

oo N / ■ \ ^ 



N=2 u.^uj)j=l 

1^1 j ^Uj 



2ti 

h(3u!c 



exp(-.|:(-l)^^^. 



■jP3 



dt 



N 



dt 



N-r ■ ■ 



t2 



dtx 



X exp < I 



N-l 

E 

3=0 



dt"[E,^{t")-E^^{t")]~^p]^,T,+, 



(F4) 



where 5j = (1) for a vertical (horizontal) jump. 
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In order to evaluate the series of integrals in Eq. (F4), we use the following technique: The upper limit t of the 
first integral is replaced by oo and for compensation the step function Q{t — ijv) is added to the integrand. Then, the 
order of integration is interchanged and the integrals are transformed to difference coordinates Tj — tj 



Like 



in the previous Sections, it is assumed that the driving frequency is large and averaging over the driving period is 
appropriate. In a final step, the 0-function in the integrand is represented as a complex integral according to 



Q{t-tN) 



1 

27ri 



dX — exp 



A t 




(F5) 



The complex integration over A can afterwards be carried out with the help of the calculus of residues for the residue 
at A = 0. After this tedious but straightforward procedure, one arrives at the final result for the averaged Markovian 
approximated rate matrix elements 



N 



-1 "J - 



^=2 {^j„-}j=l 



N 



27r 



exp 



dr Jo [p ~ sin 



exp 



(F6) 



where Jo{x) is the zeroth Bessel function. 

It is not possible to treat this complicated expression analytically. If one has to determine explicit results, the help 
of the computer is needed to evaluate the sum over all configurations {njVj}. Then, in a first step, all paths belonging 
to a fixed order N are created numerically by means of recursive programming. In the next step the sum over all the 
occurring paths has to be evaluated before one has to go to the next order by increasing N . Finally, convergence with 
respect to N has to be obtained. 

We summarize this Appendix by concluding that also the effect of the flow to weak damping can be treated by 
real-time path integrals although the expressions become much more involved. Eq. ( p^ ) constitutes the starting point 
for the study of the cross-over to the classical regime. One has to be aware, however, that the driven problem is far 
from being trivial since even a chaotic dynamics may occur. 
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FIG. 1. The states of the reduced density matrix of an M-level system. Shown are two examples of paths that travel 

between two diagonal states of the density matrix (see text). One path (solid line) connects the diagonal states O and the 
other (dashed line) travels between the diagonal states The intermittently visited off-diagonl states are marked by •. 
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a.) Representation with localized states 



b.) Discrete variable representation DVR 
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q q^, q„, q %_ % 

FIG. 2. a.) The first four localized states {q\Li) , {q\R2) of tlie static symmetric double-well potential, Eq. (^) with 
£ — s = 0, are shown in position space. They are defined according to Eq. The barrier height is chosen to be 

Eb = AU/huio = 1.4 (we use here and in the following figure captions dimensionless quantities according to the standard 
scaling defined in Eq. (A.1) ). The energy eigenvalues £i, ...,£4 axe marked by thin solid horizontal lines, b.) The corresponding 
four DVR-states are shown, i.e., {q\ai) (solid line), (g|a2) (dashed line), (gift) (dashed-dotted line), and {q\Pi} (long-dashed 
line). On the g-axis, the exact eigenvalues q^ are marked by crosses. 
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FIG. 3. The probability Picft of finding the particle in the left well as a function of time for the symmetric (e — 0) and the 
asymmetric (e = 0.08) case. Considered is the system of two doublets, i.e., M = 4. We start from an initially fully localized 
state in the left well. The barrier height is set to Eb = 1.4. The temperature is T = 0.2, the damping constant is 7 = 0.1 and 
the cut-off frequency is ujc = 10.0. For this set of parameters, the dynamics is fully incoherent and the Markov approximation 
to the GME (p8|) is rather satisfactory. 



a.) 



With driving: s=1.0, ^=0.815 



b.) 



Witli driving: s=1 .0, Q.=0.8^5 
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FIG. 4. Same as Fig. ^ in presence of resonant driving s{t) = s sm{D.t) (s = 1.0, = ujo = 0.815) for T = 0.1 (a.) and for 
T = 0.2 (b.). Insets: enlarged parts of the figures. 
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No driving: s=0, symmetric case: e=0 
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FIG. 5. Quantum relaxation rate F for the static symmetric double-well potential with barrier height Eb ~ 1.4 in dependence 
of the number M of energy eigenstates for the damping constants 7 = 0.1(*), 7 = 0.5(n) and 7 — 1.0(A). The temperature is 
chosen to be T = 0.1 and the cut-off frequency is Uc = 10.0. 
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No driving: s=0. 
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FIG. 6. a.) Spectrum of the unperturbed static (s — 0) system Hamiltonian (|l]) for a barrier height of Eb = 1.4 as a 
function of the static bias e. b.) Quantum relaxation rate F according to Eq. ( [75[ ) as a function of the static bias e for different 
temperatures T. The barrier height is Eb = 1.4 and the number of energy levels is M = 4. The bath parameters are 7 = 0.1 
and ujc ~ 10.0. 
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No driving: s=0 




No driving: s=0 
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FIG. 8. Asymptotic population P{^^ of the left well in absence of time-dependent driving as a function of the asymmetry e 
for two different temperatures (solid line) for Ad = 8. The parameters are: Eb = 1.4, 7 — 0.1 and ujc ~ 10.0. The dashed line 
marks the results obtained from a Boltzmann equilibrium distribution for the same parameters. Inset: P{^^ vs. the number M 
of energy levels for a fixed asymmetry e = 0.08 (solid line: gNICA in second order; dashed line: with Boltzmann distribution). 
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With off-resonant driving: s=0.1 , £1=0.2 




FIG. 9. Averaged quantum relaxation rate F^^ as a function of the static bias e in presence of an off-resonant driv- 
ing with s = 0.1, Q — 0.2. Shown are the results for three different numbers M of energy levels. The parameters are 
Eb ~ 1.4, T = 0.1, 7 = 0.1 and ujc — 10.0. Inset: The corresponding asymptotic populations P{^^ of the left well as a function 

of £. 
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FIG. 10. Averaged quantum relaxation rate F^^ as a function of the static bias e for three different temperatures T = 0.1 
(solid line), T = 0.3 (dashed line) and T — 0.5 (dotted line) for the resonantly driven double-doublet system M = 4. Shown 
are the cases of weak driving s = 0.1 (for T — 0.1 only) and strong driving (s = 1.0). For the remaining parameters, cf. Fig. H. 
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With resonant driving: s=0.1 , £1=0.81 5 
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FIG. 11. r''^ as a function of the static bias e for four different combinations of the number of levels, i.e., M = 6 (•) and 
M — 8 (□), and different temperatures, i.e., T = 0.1 (solid line) and T = 0.15 (dashed hne). For the remaining parameters, 
see Fig. 



With resonant driving: £1=0.815 
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FIG. 12. Asymptotic population P{^^ of the left well as a function of the asymmetry e for the temperatures T = 0.1 (solid 
line), T — 0.3 (dotted line) and T — 0.5 (dashed hne) for the double-doublet system M = 4. Shown are results for strong 
driving s = 1.0 and for weak driving s = 0.1 (for T = 0.1 only). For the remaining parameters, see Fig. ^. Inset: P{^^ vs. the 
static bias e for a fixed temperature T = 0.1 for M = 4 (solid line), M = 6 (□) and M = 8 (•). 
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With driving: il=0.815, e=0, IV1=4 




S 

FIG. 13. Averaged quantum relaxation rate F^^ as a function of the driving amplitude s for three different temperatures 
T — 0.1 (sohd line), T = 0.3 (dashed line) and T — 0.5 (dotted hne) for the driven symmetric (e — 0) double-doublet system 
M = 4. The static barrier height is Eb = 1.4 and the driving frequency is Q — luq = 0.815. The bath parameters are 7 = 0.1 
and ujc = 10.0. The asterisks * mark the findings of an exponential fit to QUAPI results (not shown). 



Willi driving: i2=0.815, e=0, T=0.1 
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FIG. 14. F^^ as a function of the driving amplitude s for an increasing number of levels. The temperature is fixed to 



T = 0.1. For the remaining parameters, see Fig, 
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With driving: s=0.1 M=4 
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FIG. 15. Averaged quantum relaxation rate F^^ as a lunction of the driving frequency Q for the driven symmetric (solid 
line) and the asymmetric (dashed line) double-doublet system M = 4. The static barrier height is Eb = 1.4 and the driving 
strength is s = 0.1. The bath parameters are T — 0.1, 7 = 0.1 and coc = 10.0. 
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Symmetric case: e=0 



b.) 



Asymmetric case: e=0.08 
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FIG. 16. F^^ as a function of the driving frequency SI for different numbers M of levels, a.): symmetric case e = 0, b. 
asymmetric case e — 0.08. The remaining parameters are as in Fig. Iisl 
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With driving: s=0.1, 8=0.08 




0.0 



0.0 



0.2 



0.4 



0.6 



0.8 



1.0 



1.2 



FIG. 17. Asymptotic population Pj^^ of the left well as a function of the driving frequency f2 for an increasing number M 
of states. The driving amplitude is s = 0.1 and the static bias is e = 0.08. The temperature is fixed at T = 0.1, the damping 
constant is chosen to be 7 = 0.1 and the cut-off is ujc = 10.0. 
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FIG. 18. (Averaged) quantum relaxation rate F*^^' as a function of temperature T for four different combinations of the 
undriven (s = 0) and the resonantly driven (s = 0.1, f2 = 0.815) case without (e = 0) and with (e = 0.08) static bias for 
the double-doublet system M = 4. The parameters are Eb = 1.4,7 = 0.1 and u!c = 10.0. Inset: Enlarged part of the low 
temperature regime for the undriven case s = 0. 
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static bias: e=0.08 
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FIG. 20. Asymptotic population P^it of the left well as a function of temperature T for an increasing number M of states. 
The static bias is e = 0.08. Shown are the case without ac-driving s = together with the results obtained from a Boltzmann 
equilibrium distribution (dashed-dotted line) and the case with off-resonant ac-driving s = 0.1,11 = 0.2. The remaining 
parameters are as in Fig. |l8|. 
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No driving: s=0 



b.) 



Willi off-resonant driving: s=0.1, Q.=0.2 
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FIG. 21. (Averaged) quantum relaxation rate r'^^^, respectively, as a function of the damping strength 7 for different 
numbers M of levels, a.): no ac-driving s = 0, b.): with off-resonant ac-driving s = 0.1, Q = 0.2. Shown are the results for the 
cases without (e = 0) and with (e = 0.08) static bias. The parameters are Eb = 1.4, T = 0.1 and ujc = 10.0. 
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FIG. 22. Example of one path consisting of three sojourns and two clusters (see text), a.) Time-resolved representation of 

the path jumping between diagonal (dashed line) and off-diagonal states, b.) The same path illustrated in the {q, g')-planc of 
the reduced density matrix. The labels qk,g'k of the reduced density matrix are not specified further. Diamonds O mark the 
visited diagonal states and filled circles • mark the visited off-diagonal states. 
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FIG. 23. a.) Exact localized states {q\Ri) and {q\R2) (full lines) of the unperturbed system Hamiltonian in position 
representation for a barrier height of Eb = 1.4. The dashed lines depict the results of the harmonic well approximation. The 
horizontal lines mark the exact eigenenergies of the double-well potential occurring always in pairs (doublets, the intradoublet 
spacing between the lower lying states is not visible on this scale). Note that the harmonic states are energetically shifted to 
the position of the exact localized states for graphical reasons, b.) Same for a barrier height Eb = 2.5. 
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